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NONSTEADY LIQUID AND GAS FLOW WITH HEAT 
ADDITION AND SHOCK PERTURBATIONS 
by Fred S. Sidransky and Margaret Marie Smith 
Lewis Research Center 

SUMMARY 

The theory of one-dimensional nonhomentropic nonsteady fluid flow is developed inde- 
pendently of a specific form of the equation of state. The method of characteristics is 
used to develop general compatibility relations which are applicable to either liquids or 
gases. By assuming specific equations of state, the classical water-hammer equations 
of Joukowski and Allievi and the nonsteady gas flow relationships of Riemann are deduced 
from the general compatibility relations. 

General numerical methods, based on the method of characteristics, are discussed 
for the two characteristic network for nonsteady liquid flow and the three characteristic 
network for nonsteady gas flow. Computer programs, utilizing these numerical proce- 
dures, are also given. 

To illustrate the versatility and to corroborate in part these methods of solving non- 
steady flow problems, three examples were selected: (1) nonsteady liquid flow, (2) non- 
steady gas flow with heat addition (or removal), (3) shock perturbations in a supersonic 
diffuser. The first example is verified by an alternate technique. The second example 
is verified in part by the Rayleigh analysis for steady flow in constant area ducts with 
heating or cooling. 


INTRODUCTION 

The analysis and study of nonsteady liquid flow and nonsteady gas flow have tended 
to diverge into two branches of fluid mechanics (refs. 1 and 2). This division has arisen 
quite naturally because of engineering requirements solely in hydraulics such as in the 
analysis of water-hammer in penstocks, and solely in gas dynamics such as in the anal- 
ysis of shock tubes and pulse jets. In the analysis of advanced systems, for example, 
rocket engines and Rankine cycle space power systems, the dynamics of the system is 



governed in large measure by the interdependent nonsteady flow characteristics of fluids 
in different phases. A representation of transients in such systems is therefore depend- 
ent on an understanding of the dynamics of fluids in both the liquid and gas phases. 

In this report, it will be shown, using the method of characteristics, that the classi- 
cal water-hammer equations of Joukowski and Allievi and the nonsteady gas flow relation- 
ships of Riemann can be deduced from a general theory for one-dimensional nonhomen- 
tropic nonsteady fluid flow. (In homentropic flow, the entropy of each fluid particle is 
equal to the entropy of any other particle in the flow region; whereas in isentropic flow, 
the entropy of each particle is constant but may be different from any other particle. 
Hence a flow may be isentropic but nonhomentropic since different particles may have 
different entropies.) 

To facilitate the use of the theory for the analysis of nonsteady nonhomentropic fluid 
flow, numerical procedures, derived from the general theory and adaptable to high-speed 
computers, are discussed. More particular attention is given to nonhomentropic non- 
steady gas flow than to liquid flow, owing to its greater complexity. (A parallel to this 
characteristic problem may be found in supersonic rotational flow. ) 

To illustrate and corroborate these numerical methods, computer programs, which 
are given in appendixes C and D, were developed at the Lewis Research Center, and the 
three examples selected are as follows: (1) nonsteady liquid flow (or water-hammer), 

(2) nonsteady gas flow with heat addition (or removal), (3) shock perturbations in a super- 
sonic diffuser. As a check on the accuracy of the numerical methods, the first example 
is verified by an alternate technique (ref. 3), and the second is checked in part by the 
Rayleigh analysis for steady flow in constant area ducts with heating and cooling. 


THEORY 

General Compatibility Relations 

The fundamental relation for continuity is given by 

v ie + p8' + J» = o (i) 

ay 3y at 

(Symbols are defined in appendix A. ) Conservation of momentum (omitting body and dis- 
sipative forces) yields 


v _3v + _3v _ 1 9P 

3y 3t p 3y 


( 2 ) 


2 


If reversible heat addition along the path of a fluid particle of fixed identity is assumed, 
the following is obtained: 


v — + — = — = _ ip 

dy at T Dt Dt 


where 


Dt at + V 3y 

(For a detailed discussion of these basic equations, see refs. 4, 5, and 6. ) The density 
p, which is some unspecified function of the entropy s and the pressure P, is expressed 
as 


P =P(s, P) 


(4) 


By equation (4), it follows that 


and 

and inasmuch as 
equation (1) becomes 

v 9P 
a 2 9 y + 

Following reference 4, let 


dp _ _ap_ _ap ( 8p as 

9y ap 9y 0s 3y 


dp _ d£_ ap ( dp 0s 

at ap at as at 



iPi® +p iz + J_^P + ^is = 0 

as ay ay 2 at as at 

cL 


(5) 


( 6 ) 


(7) 


( 8 ) 


3 



k =!(y,t) 


0) 


v = q(y, t) 


( 10 ) 


Equations (9) and (10) are used to transform equations (2), (8), and (3) to the following: 


/yJL+JL) = -2El*L-*z (v^L+m 

34 P 3 y 3| \ 3y at/ 3rj p 3y dq \ dy dt) 


( 11 ) 


^W v ii + iL\ + iz ii + isj£ / ag. + ag\ = .iP_L f v Ja+J2) 


3^ 2 \ 3y 3t/ 34 3y 34 3s \ 3y 3t, 

a 


dq 2 \ 3y 3ty 


_ Jv p _3?? _ Js / 3 t7 + _3^\ 3£ (12) 

3?7 3y 3?7 \ 3y 3t / 3s 


MLk + M\ = .MLM + ^) + ^ 

34 V 3y at/ 3*7 \ 3y 3t/ 


(13) 


From these equations, the compatibility relations and their corresponding characteristic 
directions are determined using the theory of characteristics as discussed in reference 4 
and in appendix A of reference 6 (vol. I). The compatibility relations will yield the prop- 
erties of the flow on a characteristic. The intersection of the characteristics indicates 
the position of a net point on the time-distance plane. The general compatibility relations, 
as derived in appendix B of this report, are given by 


6+P + 5!v = . ^ 


pa 5t 


5t 


3 (In p) 
3s 


(14) 


for the characteristic having a positive slope, 


1 6 " p S~v - al/ , 9 ( ln P) 


pa 5t 6t 


3s 


(15) 


for the characteristic having a negative slope in subsonic flow, and 


Ds = l /, 

Dt 


(16) 


4 



6 * 

for the particle path. The directional derivatives — and 

6t 

compatibility equations (eqs. (14) and (15)) are defined by 


r *~ 

— appearing in the general 
6t 


= ± + ( v + a) -L 

(17) 

6t at 3y 


C = A + (v -a)A 

(18) 

6t at 3y 



The characteristic slopes corresponding to the general compatibility equations (eqs. (14), 
(15), and (16)) are defined, respectively, by 


*-v + a 

(19) 

dt 


^ = v- a 

(20) 

dt 


> 

II 

£1 

(21) 


dt 


Because the compatibility equations as expressed in equations (14), (15), and (16) are not 
functions of any particular equation of state, they may be used for the one-dimensional 
nonsteady nonhomentropic flow analysis of either a liquid or a gas. Indeed, it can be 
shown that both the classical water-hammer equations of Joukowski and Allievi (ref. 2) 
and the nonsteady gas flow relations of Riemann (ref. 5) can be deduced from these com- 
patibility relations. 


Liquid Dynamics 

If it is assumed that in a liquid the percent change of density with entropy change is 
negligibly small, then it may be assumed that 

d ( l rL _g) = 0 (22) 

3s 


5 



Thus the compatibility equations become 


pa 6t 6t 


(23) 


and 


1 6~P 5~v_ 0 

pa 6t at 


(24) 


Since there are only two equations with two unknowns (viz. , the pressure P and the ve- 
locity v), the third compatibility relation (eq. (16)) may be omitted. 

The compatibility relations (eqs. (23) and (24)) may be presented in a more familiar 
form by defining the head H as the pressure divided by the specific weight of the fluid on 
the Earth's surface (here, the assumption is that the datum or reference pressure is 
zero); then 


- = g H (25) 

P C 

where g = 32. 2 feet per second squared. If the volume flow is obtained from 

q = Fv 

then the compatibility relations become 

5 + H a 6 + q 

St g c F at 

5~H _ a 5~q 

St g F 6t 
°c 

which are the Joukowski water-hammer relations in terms of head and volume flow. 

The fundamental or canonical water-hammer relations may be derived from equa- 
tions (27) and (28) . If the velocity v is assumed to be negligible as compared to the 
acoustic velocity, equations (27) and (28) become, respectively, 


(26) 

(27) 

(28) 
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( 29 ) 


_ 3H a 2 9q Q _ 9H a 9q n 
g — + 3 + ag — + * = 0 

c 3t F 3y c 3y F 3t 

9H + a2 9q_ JH . a9q =0 
c 3t F 3y c 3y F 3t 


(30) 


after expanding according to the definitions of the directional derivatives. Equations (29) 
and (30) are obviously true if 


9H _ _a^ 9q 
3t F 3y 

(31) 

9H _ __1 _9q 
3y F 3t 

(32) 


which are known as the canonical water-hammer equations (ref. 2). If the acoustic veloc- 
tiy is a constant, these canonical equations may be shown to be but another form of the 
classical wave equation 


9 2 H 

3y 2 


1 9 2 H 

a 2 3t 2 


(33) 


Thus for hydraulic systems in which the flow velocity is small as compared to the acoustic 
velocity and in which the acoustic velocity does not vary along a characteristic, the fore- 
going system of equations is applicable. This, it should be noted, does not exclude the 
possibility that different characteristics in the flow region may have different acoustic 
velocities as may exist in a duct with discrete cross-sectional area discontinuities. In a 
hydraulic system with a gradually varying area conduit and in which body and dissipative 
forces are significant, the compatibility relations may be rederived from the fundamental 
equations (viz. , eqs. (1), (2), and (3)) and put in a more general form. The canonical 
equations (eqs. (31) and (32)) are not applicable in this case 


Gas Dynamics 


ate 


The assumption of equation (22) is, of course, untenable in a perfect gas. 

dfl 11 P) . f or a g as> first consider that 
3s 


To evalu- 


7 


( 34 ) 


ds = C d(ln T) - -5- d(ln P) 

P S C J 


from thermodynamics and the perfect gas law. Inasmuch as 


C„ = 


— — 
\Y - V g c J 


(35) 


and if y is assumed constant, equation (34) may be transformed to 


ds = d(ln a) - d(ln P) 
R y - 1 


(36) 


since 


= VyRT 


(37) 


Now from equation (37) and the equation of state 


P = ^ a 2 


(38) 


Hence, instead of equation (36) 


— ds = — d(ln a) - 2 d(ln a) - d(ln p) 
R y - 1 


(39) 


is obtained. Equation (39) may be simplified to 


d( ln p) = /_l_W_!^ ds 
\Y ~ V 2 R 


(40) 


By equations (35) and (37), and since for a perfect gas 


h = C p T 


(41) 
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equation (40) may be rewritten as 


, v S C J 1 g c J 

d(ln p) = dh ds 

yR T R 


3 (In p) 

. g C J 1 | 


. _ g c J 

3s 

p ^ RT 

W 

p R 


By Maxwell’s thermodynamic relationship 

= T 
P 

it can be seen that 

d (In p) _ g c J /l - y\ _ 1 

3s R \ y / C p 

Thus for a perfect gas, the compatibility relations (eqs. (14) and (15)) become 



J_ + = *£ (y - l\ 

pa 8t 6t R \ y / 



for a characteristic with a positive slope (cf. eq. (19)) and 

- 5 ~ v - agc<T (y - A xp 

pa 6t 5t R \ y / C 

P 

for a characteristic with a negative slope in subsonic flow (cf. eq. (20)) and 

Ds = )/ / 

Dt 


(42) 


(43) 


(44) 


(45) 


(46) 


(47) 


( 16 ) 


for the particle path. 
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In the special case of homentropic flow where the entropy is not a function of the con- 

2 2 

duit location y or the time t, the Riemann variables a + v and a - v are 

y - 1 y - 1 

constants; this may be shown from the aforementioned compatibility relations by consid- 
ering the following form of equation (46): 


1 P 6 + P | 6+v _ ag c J ly - l\ ^ 
pa P 5t 6t R \ y / 

With the help of equations (36) and (38), equation (48) becomes 


17 2y \ S + (ln a) g c J S + s 

} 6 + v _ ag c J ly - l\ 

1 

4-> 

O 

i 

to 

1 

St R \ y ) 


(48) 


(49) 


Ds 

Further manipulation and remembering that \p = — yields 

Dt 


— / 2 a + v \ - ag c J ly - l\ Ds + a 6 + s 


St \y - 1 


R 


y / Dt y R St 


(50) 


For homentropic flow, the right side of equation (50) is zero. Hence, 


— (■ — - — a + v ) =0 
6t \y - 1 / 


and 


a + v = const (51) 

r - i 

along the characteristic direction given by equation (19). A corresponding derivation may 
be made from equation (47) yielding 


a - v = const (52) 

y - i 

along the characteristic direction given by equation (20). Thus, for homentropic flow, 
the compatibility relations as presented in equations (46) and (47) give the classical 
Riemann variables. 

10 
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As in reference 5, the compatibility relations for nonhomentropic gas flow may be 
made nondimensional by specifying a reference acoustic velocity a Q , which may be se- 
lected from steady- state conditions and a specified reference length y . From these, a 
reference time t may be deduced from 


‘o=? 

a o 


(53) 


Furthermore, two nondimensional parameters may be defined as 


C--Z- 


(54a) 


r = — 


(54b) 


Multiplying the compatibility relation (eq. (50)) by y Q / a Q results in 


5+/ 2 * + + 


6 t \y - 1 


Dt 


6t 


(55) 


where the nondimensional entropy S is 


g r J s 

s =_£— 

y R 


the nondimensional acoustic velocity is 



and the nondimensional velocity is 


« = — 


(56) 


(57) 


( 58 ) 


4 

is. 
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Time 


It may also be shown that the other compatibility relations in nondimensional form become 


6r \y - 1 / Dt 6t 

and for the particle path 

PS _ g c Jt o ^ 

Dt yR 


(59) 


(60) 


NUMERICAL PROCEDURES 

Each net point on the distance-time plane is determined in general by the compatibil- 
ity relations and their corresponding directions. The properties of the flow such as pres- 
sure and entropy can be found from the compatibility relations. The characteristic direc- 
tions serve to locate the position of the net point on the distance-time plane. 

Two options are before the computer in the construction of a characteristic network - 
either a free or a fixed characteristic network may be developed. In a free characteristic 
network, the location of a new or unknown characteristic net point, as net point 3’ in fig- 
ure 1, is determined from known end points A and C on the respective left running and 
right running characteristics. Because it is less laborious, the free characteristic net- 
work has been commonly utilized in hand computation (ref. 5, pp. 43 to 45). In a fixed 
characteristic network, the position of the net point on the distance-time plane is pre- 
scribed as net point 3 in figure 1. Consequently, the characteristics passing through this 

point must be determined, which means that the 
end points of the characteristics such as 1 and 2 
in figure 1 must be found by interpolation. Fur- 
thermore, in a fixed network procedure either 
the time or distance coordinates of the end points 
must be given, the other coordinate being deter- 
mined analytically. 

In this report a fixed network was utilized 
because it satisfied typical engineering require- 
ments for information as to dynamic flow condi- 
tions at prescribed locations and uniform time 
increments. The interpolation that is necessary 
for finding end points 1 and 2 as well as 4 pre- 


Particle path 

Left running // $ , 

direction (P or / I A/\ 

C, wave)^ / !/t \ \ C h wave) 

Vv/ // \ \ 

sy / I \ \ 



Right running 
direction or 


Free characteristic net 

Fixed characteristic net 


Distance — »- 

Figure 1. - Basic net point 
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sents no special difficulty with the availability of a high-speed computer. In the examples 
of this report, the time coordinate of the end points was preselected for nonsteady gas 
flow while, for the sake of convenience, the end points were prescribed at fixed distance 
coordinates for nonsteady liquid flow. 


Numerical Methods for Liquids 

Ba sic net point . - Integrating equation (27) along the left running characteristic re- 
sults in 


H + — q = C (61) 

g F a 
6 c 

where C a is the constant of integration or water-hammer variable for the left running 
wave, which is analogous to the nonsteady gas flow equation (eq. (51)). By equation (61) 
the following may be written along the left running characteristic: 


a. 


H i + 


Sc F l 


qj = h 3 + 


g C F l 


«3 


= C 


(62) 


Integrating equation (28) along the right running direction gives 


H - 


g f 
& c 


q = C, 


(63) 


where is a constant of integration or water-hammer variable for the right running 
wave, which is similar to equation (52). Furthermore, equation (63) yields 


H 


2 





(64) 


The corresponding characteristic direction for the left running wave is — = a, and for 

dv ^ 

the right running wave, — = -a, since the velocity v is assumed to be negligibly small 

dt 

compared to the acoustic velocity a. Obviously the volume flow q^ at net point 3 is 
determined by the simultaneous solution of equations (62) and (64), which is 
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( 65 ) 


Sc P l + Sc F 2 

The volume flow at net point 3 can then be substituted into either equation (62) or (64) to 
compute the head at that net point. 

Inasmuch as in liquid dynamics the location of end points 1 and 2 is preselected (as 
mentioned previously), the corresponding time coordinates of the net points must be 


t 


1 



*3-*! 

a l 


( 66 ) 


*2 



y2-y 3 

a 2 


(67) 


the distance coordinates yg and y^ being fixed. Since the acoustic velocities a^ and 
are constants, the times tj and t ^ need be computed but once. 

Summing up, the computation of Hg and q 3 may proceed as follows. From known 
values of t ^ and tg and the values of head and volume flow at those times, the water- 
hammer variables C and C. may be computed as in equations (62) and (64). The vol- 
ume flow is then defined from equation (65), and the head Hg may be computed from 
equation (62) or (64). This process is continued proceeding from one net point to another, 
net point 3 being always the net point whose position is known on the distance-time plane 
but whose head and flow are unknown. 

Since the characteristic curves may be regarded as patching curves (cf. ref. 6, 
p. 601), the patching of flow regions which are analytically different is permitted. This 
allows for the introduction of a valve, for example, at any given distance coordinate, 
which may be characterized by 


«3 


= K 


(H 3£ - H 


3r' 


H 3£ - H 3r 


1/2 


( 68 ) 


where is the valve orifice coefficient and Hg^ and Hg r are the heads immediately 
ahead and behind the valve. In this instance, equations (62) and (64) become 
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(69) 


H i + 


a * a l 

— Ho/) ■+ * Q.g — C, 


*C F 1 


q l - H 3£ 


9»X 


c J l 


Ho 




c 2 


q 2 - H 3r 


g F 


q 3 = C t 


c 2 


(70) 


Equations (68), (69), and (70) may be solved simultaneously, and using the Quadratic 
Formula (ref. 7) results in 


q 3 


= K 


2 Ok - V < c a - 


c - c, 

a b 


-1 + 


{ 


4 C - CL 
a b 


Kok - V 2 


(71) 


where 



and 




(72) 


(73) 


This should be compared to reference 3, page 10 in which the field method is used instead 
of the lattice point method of this report. These two methods of solving the partial differ- 
ential equations of nonsteady (and supersonic) flow are discussed in reference 6, 
pages 491 and 492. 

Left and right boundaries . - For the left boundary, only the right running compatibil- 
ity relation (eq. (64)) need be used. If the end point of the right running wave is repre- 
sented by net point 2, appropriate values of head and flow, Hg and qg at time tg, 
define the water-hammer variable C^. If net point 3 is chosen as the left boundary point, 
there are two unknowns (viz. , the head and flow at the boundary Hg and q^). Hence, 
either the head or the flow at the boundary must be given; for example, the head at the 
outlet of a reservoir is typically assumed to be constant. Of course, the boundary may 
be a valve in which case the functions describing valve performance together with the com- 
patibility relation can be solved either analytically, if possible, or numerically. 
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For the right boundary, only compatibility relation (eq. (62)) is required, and the 
past reference time is clearly tj. For some assumed boundary condition, the procedures 
are quite similar to the left boundary discussed previously. 


Numerical Methods for Gases 

Basic net point . - In nonsteady nonhomentropic flow, the procedure for the basic net 
point (cf. fig. 1, p. 12) is considerably more difficult because of essentially two compli- 
cations: (1) the particle path or third characteristic cannot be neglected and (2) the char- 
acteristic slopes are not constants. In the first instance, the particle path cannot be ne- 
glected since the compatibility equations for a perfect gas (viz. , eqs. (55) and (59)) depend 

DS 

on the value of the co-moving or substantial derivative — , a fundamental parameter for 

Dr 

the third characteristic. Thus to find the flow conditions at net point 3 in figure 1, three 
compatibility equations instead of two must be solved. (Net point 3 always identifies the 
"later” net point whose position is known on the distance-time plane but whose flow prop- 
erties are unknown.) Secondly, owing to the compressibility of a perfect gas, the acoustic 
velocity cannot be assumed to be a constant. In proceeding from net point 1 to net point 3, 
for example, the acoustic as well as the flow velocities at each net point may be different, 
requiring an approximation of the characteristic slope (cf. eq. (19)). Although the posi- 
tion of net point 3 on the distance-time plane is preselected, it does not follow that the 
location of net point 1 on the base line is immediately evident as is the case in liquid dy- 
namics. Rather, due to the variability of the acoustic and flow velocities, the location of 
net point 1 becomes part of the problem. 

The difficulties encountered and the solution to the basic net point procedure in a 
fixed network may be more readily understood by considering the following example. In 
figure 1, at points A, B, and C on the initial line all parameters are assumed to be known. 
Assume that the end points 1 and 4 are between A and B, and that end point 2 is between 
B and C. Also assume that reasonable guesses for the nondimensional acoustic velocity 
at net point 3 and the nondimensional velocity are st-^ and respectively, 
which are the flow velocity and the acoustic velocity at B. By linearly interpolating for 
j/j and between A and B, the correct location of net point 1 on the base line ABC 
will be that point whose characteristic slope will pass through the preselected net point 3 
and net point 1. The characteristic slope in this instance is given by 


t 3 ~ T 1 2 

? 3 -?l (% + %) + (j*i + sty 


(74) 
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which is an inverted finite difference form of equation (19). It may be shown by the appli- 
cation of the elementary principles of analytical geometry that the location of net point 1 
£ j may be expressed in general by the quadratic 

a jfil + Pjfil + ( Pl = 0 ( 75 > 


where 


and 


a l ~ C alf C uU 


&$_ = C al! C n2£ “ C nM C ul£ ' 2 


cp fL -2t >z - C nU C n2£ 


'al£ 


T A " T B 
? A _ ? B 


"ulf 


1 + C 


a If 


(? 


A 


;„> 2 + < t a - t b > 2 


1/2 


[<«B - *A> + <^B - 4>] 


C nlt ' r 3 - t A + C alt ! A 


C n2J! ” *3 + 4 + 4 + 4" C alt ? A 


Likewise, for the location of netpoint 2 £ 2 


«r?2 +/3 r?2 + < ^r = 0 


and for the slope (cf. eq. (20)) 


(76a) 

(76b) 

(76c) 

(77 a) 

(77b) 

(77 c) 
(77d) 

(78a) 
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2 


(78b) 


where 


t 3 “ t 2 


" ^2 % + % ~ **2 ~ ^3 


and 


a r ~ ^alr^ulr 


^r _ ^alr^n2r _ ^nlr^ulr " ^ 


" 2? 3 " C nlr C n2r 


(79a) 

(79b) 

(79c) 


'air 


t B ' T C 
? B " S C 


(80a) 


■'u lr 


1 + C 


air 


^B " + ^ t B " T C ) _ 


1/2 


[<*C ” " ^B>] 


(80b) 


'nlr = T 3- T C + C alr^B 


C n2r " *3 ” ^3 + - *^B " C ulr ? B 


(80c) 

(80d) 


For the location of £ ^ the end point of the particle path on the initial line, a similar 
quadratic is the result 


a t. + 0 t.+(p =0 

m s 4 F m s 4 


(81a) 


and for the slope (cf. eq. (21)) 


t 3 ~ t 4 _ 2 

? 3 _? 4 % + 4f 4 


(81b) 
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where 


01 m ^alm^ulm 


^n2m^alm " ^ulm^nlm " 2 


" ^3 " C nlm C n 2 m 


(81c) 
(8 Id) 
(81e) 


and 


'aim 


t A ~ T B 
? A ' 


(82a) 


'ulm 


1 + C 


aim 


" ^ a } + ^ t B ” r A^ 


1/2 


(♦n - «a) 


(82b) 


'nlm r 3 r A + ^alm^A 


'n 2 m = *3 + *A “ <W A 


(82c) 

(82d) 


These quadratics (viz. , eqs. (75), (78a), and (81a)) of course, may be solved for 
and £ 4 , respectively, by the well-known Quadratic Formula where the correct sign is 
given by the solution for which is to the left of and nearest to £ 3 , and for £ 2 which 
is to the right of and nearest to £ 3 , and for £ 4 which is to the left of and nearest to £3 
if (< 8(3 + £ & 4 )/2 > 0. Further it will be seen that, if the slope of the base line is zero, then 


a H = a r = a m = 0 


and the solutions to equations (75), (78a), and (81a) may be represented by 


Cl = 


n 


(83a) 
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(83b) 


? 2 



54 



(83c) 


respectively. 

If the estimated location and associated flow parameters of the end points of the char- 
acteristics are determined from the assumed values of si g and ‘fcg, a better approxima- 
tion of the values of j/g and ^ can be made by the Method of Iteration for Simultaneous 
Equations (ref. 8). A new value of may be determined from the assumed values of 
j/g and If the compatibility relation (eq. (60)) is expressed in finite difference form, 
the result is as in reference 5 


S 3 =S 4 + i? (T 3- T 4) 


(84) 


The time increment Tg - is obviously known since the basic network is for fixed time 

and distance coordinates. The nondimensional entropy S A is found by linearly interpolat- 

’ DS 

ing along the base line ABC at location £ If the rate of entropy increase — is 

* Dr 

assumed to be known, then Sg is defined. With the interpolated values of the flow param- 
eters at locations £ j and the Riemann variables at location £g may be found by the 
finite-difference form of equations (55) and (59), namely, 


■*. 3 + 


+ ‘i'l + (y - !) y i3 ^ ( T g - T j) + J^ 13 (Sg - Sj) 


(85a) 


and 


si n - 


r - i 


*2 " ^2 + ^ ^23 ^ t 3 ‘ T 2 ) + ^23 ^ S 3 ” S 2 ^ 


(85b) 


where the double subscript represents an average value (e. g. , j/j g = (si^ + ^/g)/2). If 


0 3~^Zl' s/ 3 + % 


(86a) 
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and 


^3 ~ ^77 y 3 ' *3 


( 86 b) 


Then ^ is determined by 


*3 " | <**3 ' •%> <«> 

a 

By following the Method of Iteration for Simultaneous Equations, the entire system 
is recalculated with the new value of 45 fg and the previous value of this includes 
finding new locations for an< ^ £ 4 * With the new interpolated values of the flow 

parameters at these locations, equations (84), (85a), and (85b) are recalculated; is 
then computed with the new value of < 2 ^ by 

■*3 - (Ss + ^ (88) 

This completes a single iteration. The process is repeated until the values of and 
4fcg converge within some desired tolerance. With a reasonable tolerance for engineer- 
ing calculations, this method has rarely required more than two iterations for the com- 
putation of any basic net point in the examples of this report. 

If, in a problem, net points 1 and 2 have different gas properties, the basic net point 
procedure is not significantly altered except that in equations (85a) and (85b) the ratio of 
specific heats y at corresponding net points must be altered to correspond with the prop- 
erties of the fluid. 

Left boundary . - If a constant pressure Pg and a constant entropy S 3 at the left 
boundary are assumed, the nondimens ional acoustic velocity at the boundary may be 
derived from equation (36) to yield 


•*3 = 


P 3\ (y -l)/2 y e [S 3 (y- 1 )/ 2 ] 

To) 


(89) 


where the nondimensional entropy at the boundary S 3 is determined from equation (56). 
The reference pressure P Q must be consistent with the reference acoustic velocity 
(cf. eq. (54)) and some assumed reference density p Q ; thus 
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as in equation (38). A solution for the nondimensional velocity at the boundary #g may 
now be found by applying the Method of Iteration (cf. ref. 8 ) . Since the location £ g and 
the time Tg at the left boundary are known, the location of on the base line may be 
found by equation (83b) if is assumed to be a good initial guess for # 3 , the nondi- 
mensional velocity at the boundary. If an initial location for ^ * s chosen, a linear 
interpolation along the base line may be made for st^, and Sg. Only a single com- 
patibility relation (i. e. , eq. (85b)) is necessary to define the Riemann variable inas- 
much as only a positive direction of flow is assumed. Then ♦g is given by 

« 3 ’—* 3-®3 < 90 > 

Y “ 1 

remembering that si g is a constant as in equation (89). If the previous value of ^ is 
different from the initial guess, a few iterations are ususally sufficient for convergence 
to some desired tolerance level. 

Later, in the discussion of the supersonic diffuser example, special procedures are 
covered if the left boundary consists of a supersonic diffuser with a normal shock down- 
stream of the inlet. 

Right boundary . - If the right boundary is choked, the ratio of the duct area F to the 
throat area F* is given by (ref. 6 , pp. 85 to 86 ) 


_F 

F* 




(y+l)/ 2 (y-l) 


(91) 


where Mg is the subsonic Mach number just upstream of the choking orifice at the right 
boundary. Hence, for a specified F/F*, Mg is constant. For positive flow, two com- 
patibility relations (eqs. (84) and (85a)) and their corresponding end point £ j and £ 4 
together with the boundary condition are necessary to define flow conditions at the bound- 
ary. If is used as an initial guess for j/g, 45fg is computed by 

* 3 = jafgMg (92) 

The locations and are defined by equations (83a) and (83c). Interpolating along 
the base line for sty and Sj at and for S 4 at £ 4 , the Riemann variable 0^ 


22 


is known by first determining Sg from equation (84) and substituting this value together 
with the aforementioned parameters in equation (85a). The acoustic velocity is then 
given by 


J* 3 = - (93) 

— — + M„ 

y- i 3 

This value is then compared to the initial guess for to assess whether additional iter- 
ations are necessary. 

Obviously the right boundary condition is not limited to a choking orifice nor is the 
left boundary limited to a constant pressure. Had the boundary conditions been reversed 
or other boundary conditions prevailed, similar methods as described previously would 
apply. 

Special net point cases . - Clearly, as in figure 2, there will be instances in the com- 
putation of the basic net point 3 when the end point of the & characteristic, net point 1, 
or the end point of the Q. characteristic, net point 2, exceed their respective boundaries. 
In those instances the boundary values must be used in place of the fictitious net points 1 
and 2. Consider the case of a net point 3 where £ j, the end point of the & character- 
istic, falls beyond the left boundary. To find the intersection of the & characteristic 
with the boundary, the location of the fictitious net point 1 is first computed by equa- 
tion (83a). The intersection of the characteristic passing through net points 1 and 3 with 



Figure 2. - Special net point cases. 
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the boundary is then found by 


T lim = r 3 " ^3 " ? lim> ~ < 94a > 

? 3 ‘ ? 1 

where T 1 irr) is the nondimensional time at the boundary and g l irn is the location of the 
boundary itself. The 21 characteristic which intersects the limit point on the left bound- 
ary as in figure 2 is determined by the methods discussed in the section on left boundary 
procedures. The values and S lim become in effect the end point flow 

parameters of the left running wave; the line connecting the limit point with net point 2 
is used as the base line on which the end point of the particle path, net point 4, is located. 
Points A and B are the limit point and net point 2, respectively (see fig. 2); quadratic 
equation (eq. (81a)) must be solved since the slope of the base line is no longer zero. 

Similar methods apply to the right boundary where, instead of equation (94a), there 
is obtained 


T lun - t 3 - «3 * film* T~T < Mb) 

? 3"^2 

Another difficulty arises when base point 1 or 2 is within the boundaries but falls to 
the left of point A or to the right of point C, respectively. In this case if net point 1 falls 

beyond A, as in figure 2, then interpolations may be performed between A 4 and B 4 . Simi- 

w 2 2 
larly, if net point 2 is beyond C, then interpolation between B and C may be necessary. 


Orangization of Network Calculations 


There are a number of ways in which a characteristic network can be organized. 

Basically each net point must be identified by 
some enumerative order; at each identified 
net point, a scheme must be devised for se- 
lecting the correct end points of the character- 
istics for that net point; and lastly, there must 
be a method of proceeding from one net point 
to another. The enumerative order selected 
in the computer program for nonsteady non- 
homentropic gas flow is shown in figure 3. 

The known base points of net point 5, 2 for ex- 
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ample, are net points 4, 1 and 5, 1 and 6, 1 (between which interpolations will probably be 
necessary) which correspond to the base points A, B, and C of figure 1 (p. 12). The first 
subscript, it will be seen, indicates the distance coordinate and the second, the time co- 
ordinate. The procedure of going from one net point to another is nothing more than pro- 
ceeding according to some numerical order along a specific time coordinate (fig. 3). 

Evidently both problems in water-hammer or gas dynamics may be solved on the 
same network using the same procedures and in fact using the same computer program if 
the general compatibility relations (eqs. (14), (15), and (16)) are employed in a region free 
of compression shocks. But inasmuch as the slopes of the characteristics are constants 
in liquid dynamics, a simpler water-hammer program may be devised (cf. appendix D 
or ref. 3). 


EXAMPLES 


Three examples have been chosen (corresponding to the three parts of fig. 4). The 
first two have been selected to corroborate the numerical methods of this report (viz. , 
an example in water-hammer or liquid dynamics which can readily be checked against the 
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(a) Hydraulic duct perturbed by periodic variation of inlet orifice. 
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(c) Supersonic diffuser. 

Figure 4. - Configurations for examples. 
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methods of reference 3 and secondly, an example in gas dynamics, the cooling of a per- 
fect gas flowing in a duct of constant cross-sectional area). A partial verification of the 
latter example is possible inasmuch as at steady- state conditions the results of the non- 
steady flow analysis should be consistent with the Rayleigh analysis for cooling or heating 
of a constant area duct (ref. 6, ch. 7). The third example, shock perturbations in a 
supersonic diffuser, has been selected to show how the analytical methods of this report 
may be applied in an approximate manner to the dynamics of a supersonic inlet, the de- 
sign of which is of current interest. 


Example 1: Hydraulic Duct Perturbed by a Periodic Variation of Inlet Orifice 

In the configuration shown in figure 4(a), three orifices are inserted at equal intervals 
of 17 feet in a constant cross-sectional area duct with an inside diameter of 7/8 inch. 

The values of the orifice coefficients are given by 

a n = %I = *rv = 2 * 456 ftl/2 /sec (95) 


The subscripts n, m, and IV denote the location of these orifices in the duct (cf. fig. 
4(a)). The orifice coefficient St is equal to the orifice coefficient divided by the 
duct area (cf. eq. (68)); hence, the velocity v at each orifice can be presented by 


v = £ 


(H £ - H r ) 
|h £ -n r \ 1/2 


(96) 


where and H r are the heads immediately ahead or upstream and behind or down- 
stream of the orifice. At the left and right boundaries are additional orifices represented 
by = 0. 65 and dt y = 0.414 (foot) ^^per second, respectively. Initially, the head just 
downstream of the first orifice is 260 feet, the head just upstream of the right boundary 
is 238.6 feet, and the velocity through the duct is 6. 5 feet per second. The acoustic 
velocity is set at 3800 feet per second throughout the duct. Perturbations are introduced 
into this simple system by varying the orifice coefficient at the left boundary accord- 
ing to the function 


*r»o* *unp Bin “* 


( 97 ) 


where & n = 0.65, = 0.5, and u> = 2nf where f = 70 cps. The time t is real 

v cLilip 
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Nondimensional time, t 


(a) Head at hydraulic line discharge plotted against nondimen- (b) Head immediately downstream of hydraulic inlet plotted 

sional time. against nondimensional time. 

Figure 5. - Nonsteady liquid flow example. 

time in seconds. 

Under these conditions the water-hammer programs of this report and of reference 3 
were used to compute the dynamics of this line. The results were practically identical, 
which may be verified by inspection of figure 5, which represents a short period or time 
slice of the entire transient. 

Example 2: Cooling Along a Constant Area Duct 

As seen in figure 4(b), a duct of constant cross-sectional area, 5 feet in length, 
choked at the right boundary, and having a constant pressure left boundary, contains, 
prior to any transients, a perfect gas flowing at the same Mach numer of 0. 384 at each 
location upstream of the choking orifice. Initial values of pressure, specific weight, and 
temperature (assumed to be constant at every upstream location) are P = 22.74 pounds 
per square inch, W = 0. 06977 pound per cubic foot, and T = 880. 57° R. A reference 
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Specific weight, W, Ib/cu ft Pressure, P, Ib/sq in. 




Figure 6. - Transients due to cooling of constant area duct Real time, t(0. 8937x10“ 2 ) second; ratio of specific heats, 
1. 4; gas constant, 53. 3 feet per °R; -j^ = -o. 5. 
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(d) Mach number transients. 

Figure 6. - Concluded. 




length y Q of 10 feet and a reference acoustic velocity a c of 1119 feet per second were 

selected. By equation (54) the reference time t becomes 0. 00894 second, and real time 

in seconds is the nondimensional time multiplied by this factor. 

A transient is initiated in this system by instantaneously imposing and sustaining a 

US 

defined rate of cooling, namely, a — = -0. 5 (cf. eq. (60)) along the duct. Periodic 

Dr 

oscillations follow the transient as can be seen in figure 6; after these oscillations have 
been spent, the flow settles about altered steady-state conditions at each location, quite 
different from the initial conditions prevailing before the transient. 

It is interesting to compare, as a check, the final steady- state conditions with the 
Rayleigh steady- state analysis (cf. ref. 6, ch. 7). After the transient has been com- 
pleted, the temperature at the left boundary is 880. 57° R and the Mach number is 0. 443. 

At the right boundary or discharge of the duct, the Mach number is 0. 384, with a corre- 
sponding temperature of 733. 79° R. With the help of table B. 5 in reference 6, pages 628 
and 629, the critical temperature ratios are (adopting the nomenclature of ref. 6) 

| = 0. 69505 (98a) 

M=0. 443 

and 

= 0. 58282 (98b) 

M=0. 384 




by interpolation. Thus the temperature ratio is 



M=0. 443 


( T *)m=0. 384 


0.69505 
0. 58282 


1. 193 


(98c) 


which compares within 1 percent of the temperature ratio of the values computed from the 
nonsteady flow analysis. A similar check may be made for pressure ratio, and it may be 
shown that excellent agreement is obtained. In addition, it may be seen that the frequency 
of the oscillations, about 59 cps, approximates the quarter wave resonant frequency for 
the duct. 
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Example 3: Perturbation of a Normal Shock in a Supersonic Diffuser 


In figure 4(c) (p. 25), a supersonic diffuser followed by a constant area duct may be 
seen. In this example the same dimensions and initial conditions prevail in the duct as 
in the previous cooling example. At the right boundary the system is perturbed by vary- 
ing the choking throat area with time (fig. 7(a)). In place of the constant pressure bound- 
ary, there is a supersonic diffuser with a diffuser exit to inlet area ratio F /F. of 

ex' in 

1.6. The Mach number at station zero (the inlet of the duct or diffuser exit) M is 
0. 384, and the Mach number at the diffuser inlet M Jn is 1. 5. In this example, the 
quasi-steady methods of reference 5, pages 99 to 103 are adopted. Hence the weight 
flow by continuity and the stagnation temperature by conservation of energy are assumed 
to be invariant from the diffuser inlet to the diffuser exit. By these assumptions if 


M 


D ex = 


ex 


1 + Z^M 2 e Y r+1)/2(r ' 1) 

2 ®*/ 


(99a) 


D in = 


M in 




in , 


(99b) 


the nondimens ional entropy increase across the shock S - S. is given by 

ex in 


s - s. 

ex in 


_ g c J ( S ex ~ S in) _ 1 ln / F ex D ex\ _ 1 ln / P T, ir 


r R 


^ F in D in/ * \ P T, 


( 100 ) 


as in reference 5 and reference 6, pages 118 to 121. For the specified M , M. 

ex in 

and F ex' F in’ S ex “ S in = °* ° 948 * The supersonic Mach number just upstream of the 

shock M may be found by solving the following equation for M : 

ss 


ex 


S. = 
in 


y - i 


ln 


+ rjLl 


(r + 4 )Mg S 


Y + 1 


ln(—^¥— 

y(r - 1) \y + l 


M* - 
ss i 

y + 1 


( 101 ) 


The ratio of shock location area to inlet area F__/F. is then known by 

bb in 
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Pressure loss, APy/Py Critical area ratio, F/F* 



(a) Periodic variations In critical area ratio in supersonic dif- (b) Change in ratio of area at shock position to inlet area plotted 



(c) Pressure loss plotted against nondimensiona! time. (d) Pressure perturbations in supersonic diffuser example at 

diffuser exit, at 3 feet, and at 5 feet 

Figure 7. - Supersonic diffuser example. 
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(102) 




where 


D ss “ 


f l +■ 


M 

ss 

.1 2 x (y+i)/2(r-i) 

■±m: ' 


ss 


(103) 


The initial ratio of shock location area to inlet area for the previously specified conditions 
is 1. 034. The initial subsonic Mach number just downstream of the shock M gu j ) may be 
found by the use of Prandtl*s equation 


M 


* _ i 

sub M * 
ss 


(104a) 


(cf. eq. (5. 17a) of ref. 6) where the dimensionless velocity M* is determined by 


*2 _ (y + 1)M^ 


M = 


(104b) 


2 + (y - 1)M" 


For this example, the value of the subsonic Mach number is 0. 654. 

To determine conditions in the diffuser, an iterative procedure is necessary. Ini- 
tially the nondimensional stagnation acoustic velocity which it will be noted is a con- 
stant, may be computed from 




= -”'ex( 1 + 


m 2 


1/2 


ex 


(105) 


using the initial values of and M. . Next, an estimate of .2 is used to determine 

ex ex ex 

the acoustic velocity at the diffuser exit .«/ at any later time from (ref. 5, p. 64) 

cx 


A 


ex 


■^ex 



(106) 
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which was derived by conservation of energy and the Zb compatibility relation. The non- 
dimensional velocity is then given by 


♦ = 

ex 


- 1 


^ex 


ex 


(107) 


The end point of the Zb characteristic and the corresponding parameters *2* ^ and 

Sg can be ascertained by methods discussed in the left boundary procedures. At this 

point the nondimensional entropy increase may be calculated by equation (100), and the 

entropy at the diffuser exit S may be deduced. The Riemann variable at the diffuser 

exit Zb may then be calculated by 
cx 

*«-^ + ^,«( s «- s a) < l08a > 

where J* 2 i ex = (-“2 + J 'ex >/2 - and 

S 2 = -~^ 2-*2 < 108b ) 

DS 

Equation (108a) is a finite difference form of equation (59) since in this example — = 0. 

Dr 

The Riemann variable Zb 2 * s defined by the interpolated values of and ‘tg 
at the end point £ 0 of the Zb characteristic. The Riemann variable computed from 
equation (108a) may be compared to the estimated value, and unless there is agreement 
between the initial guess and final calculation of Zb within a reasonable tolerance, ad- 

GX 

ditional guesses will be necessary. After these iterations have been completed, the shock 

Mach number M__ and the shock location represented by F__/F. may be computed 
ss ss in 

from equations (101) and (102), respectively. Moreover, the pressure loss AP^/P^ is 
known from 



F. D. 
in in 


F D 
ex ex 



(109) 


At the right boundary, the system, it will be remembered, is perturbed by varying 
the choking throat area with time. Clearly as the duct to throat area ratio at the right 
boundary decreases (i. e. , the throat area is being enlarged) the Mach number at the right 
boundary increases or an expansion wave is initiated at this boundary. When this wave 
finally arrives at the left boundary, the shock moves forward as seen in figure 7(b) in- 
creasing the supersonic Mach number ahead of the shock and the pressure losses (see 
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fig. 7(c)). As the choking orifice at the right boundary returns to its former position, 
compression waves resulting from this motion cause the shock to return nearly to its 
former position. Had the compression waves been too strong, the shock would have been 
blown out of the diffuser altogether. Pressure oscillations at specified locations along the 
duct due to this motion are shown in figure 7(d). 


CONCLUDING REMARKS 

A theory, numerical methods, and computer programs for one-dimensional, non- 
steady, nonhomentropic fluid flow have been presented; and it has been demonstrated that 
from the same general compatibility relations, nonsteady liquid flow and nonsteady gas 
flow with heat addition and shock perturbations may be analyzed. Moreover, in two of the 
examples selected, verification of the computations by alternate methods has been shown. 

The theory, presented herein, does not preclude the possibility of analyzing a fluid 
in which the properties (e. g. , molecular weight) of the fluid on each side of a defined 
interface are different. Such an instance may arise in a combustion chamber due to the 
initial injection of fuel, or during the expulsion of air from an open-ended duct as, say, 
nitrogen gas issues into the duct from a high pressure bottle at the opposite end. 

In the formulation of the general compatibility relations, body and dissipative forces 
were omitted together with the variation of cross-sectional area in a conduit. The funda- 
mental equations (i. e. , conservation of energy, momentum, and continuity) may be ex- 
tended to include these effects, and new compatibility relations may be derived which will 
have a broader application in nonsteady liquid and gas flow problems. Further, since no 
equation of state has been specified, the theory is not limited only to the dynamic analysis 
of liquids and gases. The state of the fluid may be described by a table of values such as 
a steam table. The theory, thus, may be applied to the dynamics of homogeneous two 
phase flow with heat addition. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 8, 1966, 

120-27-04-27-22. 
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APPENDIX A 


SYMBOLS 


j/ 

nondimensional acoustic velocity 

q 

volume flow, cu ft/sec 


defined in eq. (57) 

R 

gas constant, ft-lb/(slug)(°R) 

a 

acoustic velocity, ft/sec 

3t 

orifice coefficient (see eq. (96)), 

C a 

water-hammer variable (see 


ft 1/2 /sec 


eq. (61)), ft 

S 

nondimensional entropy defined in 

C b 

water-hammer variable (see 


eq. (56) 


eq. (63)), ft 

^a 

water-hammer parameter (see 

S 

specific heat at constant pres- 


eq. (72)), sec/sq ft 


sure, Btu/(lb)(°R) 

•*b 

water-hammer parameter (see 

D 

parameter (see eqs. (99a) and 


eq. (73)), sec/sq ft 


(99b)) 

s 

specific entropy, Btu/(lb)(°R) 

F 

cross-sectional area, sq ft 

T 

temperature, °R 

f 

frequency, cps 

t 

time, sec 

Sc 

acceleration due to Earth's 

V 

nondimensional velocity defined in 

H 

h 

J 

gravitational field, 32. 2 
ft/sec 2 

head, ft 

enthalpy, Btu/Lb 

mechanical equivalent of heat, 
778. 26 ft-lb/Btu 

V 

W 

y 

y 

eq. (58) 

velocity, ft/sec 
specific weight, lb/cu ft 
conduit length, ft 
ratio of specific heats 

K v 

orifice coefficient (see eq. (68)), 

f 

nondimensional distance defined 

M 

ft 5 / 2 /sec 
Mach number 

V 

in eq. (54a) 

parameter (see eq. (10)) 

ss 

shock Mach number 

a 

parameter (see eq. (9)) 

p 

? 

pressure, lb/sq ft or lb/sq in. 
Riemann variable (left running) 

p 

T 

density, slug/cu ft 

nondimensional time defined in 
eq. (54b) 

Q 

heat per pound of fluid, Btu/Lb 

* 

parameter (see eq. (3)) 

2 

Riemann variable (right running) 

0) 

frequency, radians/sec 
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Subscripts: 

A base point (see fig. 1) 

a refers to left running wave in 

liquid dynamics (see fig. 1) 

amp amplitude (see eq. (97)) 

B base point (see fig. 1) 

b refers to right running wave in 

liquid dynamics (see fig. 1) 

C base point (see fig. 1) 

ex diffuser exit 

in diffuser inlet 

lim limit 

£ left or upstream of fluid par- 

ticle 

m refers to fluid particle 

o reference 

r right or downstream of fluid 

particle 

ss supersonic 


sub subsonic 

T stagnation 

0 see eq. (97) 

1 end point of left running wave 

(see fig. 1) 

2 end point of right running wave 

(see fig. 1) 

3 intersection of left running wave, 

right running wave, and par- 
ticle path (see fig. 1) 

4 end point of particle path 

l, n, see fig. 4(a) 

m, iv, 
v 

Superscripts: 

1 alternate base point (see fig. 2) 

2 alternate base point (see fig. 2) 

* refers to free characteristic net 

point (see fig. 1) 

* signifies state at which M = 1 

(see ref. 6) 
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APPENDIX B 


DERIVATION OF GENERAL COMPATIBILITY RELATIONS 


ap 2 c 

First, equations (11), (12), and (13) must be solved for ~ and — according to the 

a v 

methods of reference 4. The other unknown is omitted since it will yield no addition- 
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al compatibility relations (cf. ref. 6, pp. 974-977). The denominator in every case is 


Hence, 


in 

p 3y 

viS + ii 
3y at 

0 

±(v*± + m 

a 2V 3y at/ 

3y 

sp / V M + ii 
3s \ 3y at. 

0 

0 

v-M + M 
3y at 

V 3y at/ 

7ag\ 2 i 
V 'a 2 

: 

V ay 3t/_ 


(Bl) 


(B2) 


One solution is 


v®i+.M = o (B3a) 

3y 3t 


Hence, 



dt 


and the second solution is 


a 2 (k) 2 = Lk+m) 2 

\3y / \ 3y at/ 


(B3b) 


(B4a) 
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(B4b) 


Hence, 


Now if 


solve for 


If II = - 
at 


± a ii = vil + ii 

3y 3y at 


^ = v + a 
dt 


^y = v -a 

dt 


(B4c) 

(B4d) 


x = -iZila _ Ix^la + la\ 

dt] p d y dt]\ dy dtj 


(B5a) 


Y = 


i£ by 

3£ 


I J Lfyia+M . 
i a 2\ ay at/ 

. |v p _|2 _ dsL 
dt] dy drj\ 

la +la\l£ 

3y at / 3s 

(B5b) 

z = 

dt] 

(r*n + la\ + * 
V ay at / 


(B5c) 

X vi! + ii 

ay at 

0 



y 

3y 

ds\ dy at/ 

= 0 

(B6) 

Z 0 

vM + li 

ay at 


- 


r JL, this determinant vanishes. If, as in reference 4, the following is assumed: 

dy 


= -(v + a) H 

at dy 


(B7) 
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which is implicit in equation (B4b), then 


X 



0 


Y -ifiaii 

3y 3s 3y 


= 0 


-all 

3y 


or 


(f - ) 2 (" paX ~ ^ + a2 = ° 

Dividing by a results in one possible solution to equation (B8), that is, 

-pX - aY + a Z = 0 
3s 


(B8) 


(B9) 


(BIO) 


along a line on the y, t plane having the slope given by equation (B4c). If the following is 
assumed: 


JL = -(v - a) (Bll) 
3t 3y 

instead of equation (BIO) 

-pX + aY - a-^ Z = 0 (B12) 

3s 

is obtained along a line on the y, t plane having the slope given by equation (B4d). By 
solving equations (11), (12), and (13) for — the compatibility relation along the particle 
path is found from 
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19 1 

3-3 x 

p 9y 

QJ 

<< 

QJ 


p-5i Y 

a 2V 8y 3t/ 

9y 


= 0 


0 


0 z 


or expanding along the third row 



since along the particle path 

vil + M = o 
ay at 


A possible solution to equation (B14) is 


Z = 0 


Therefore, equation (B5c) results in 


M(vin+M] = ^ 
^ \ ay at / 


(B13) 


(B14) 


(B15) 


(B16) 


(B17) 


but 


v ^z + ia = 2z? (bis) 

ay at Dt 

or 

i® 5Z? = 5^ = xjy (B19) 

drj Dt Dt 

If there is no heat addition to the system, then 
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5£ = o 

Dt 


(B20) 


or 


s = const 


(B21) 


as in equation (15) of reference 4. 

At this point, the left and right running compatibility equations of reference 4 are ex- 
tended to accommodate nonhomentropic flow in which there may be not only shock waves 
but also heat addition. From equations (B5a), (B5b), and (B5c), the left running compat- 
ibility relation (eq. (BIO)) may be presented in the following form: 

1 JP/iia + via + ia\ + p iz&ia + via + ia \ 

a dt] V 3y 9y 9t / dq \ 9y 9y 3t / 

+a ^^ v ^ + a 9s 9^ 9p. a 9p 9s v 9rj _ a 9p 9s 92 + aj// 9p = 0 (B22) 

9?7 3s 9y 9»7 9t 9s 9s dt] dy ds dq 9t 9s 


Now if 


and 


6 + 9 / x 9 

— = — + (v + a) — 

6t 9t 9y 


6“ 9 / . 9 

— = — + (v - a) — 

6t 9t 9y 

then equation (B22) may be expressed by 

1 9P5 + 77 + 9v6 + 77__^, 9(lnp) 

pa 9/7 6t dt] 6t 9s 


or 


1 6+P ., 6 + y__ flV , 9(lnp) 
pa 6t 6t 9s 


(B23a) 


(B23b) 


(B24) 


(B25) 
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which is the compatibility equation for a fluid (gas or liquid) including heat addition in the 
left running direction (cf. fig. 1). A similar development may be made for the right run- 
ning direction. Such a development yields 

±£P_ (Tv = _ a ^ 3(ln_p) (B26) 

pa 6t 6t 3s 

which is the compatibility equation for a fluid (gas or liquid) including heat addition in the 
right running direction (cf. fig. 1). 
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APPENDIX C 


COMPUTER PROGRAM FOR NONHOMENTROPIC NONSTEADY GAS FLOW 

In the following, the input or required data will be presented for this program so that 
it may be readily employed. The input of this computer program contains several options 
not utilized in either the example for nonsteady gas flow with cooling or the supersonic 
diffuser example. This will be made clear in the following outline of the input required 
for this program. 


Input Variables and Explanations 

First input card (corresponding to read statement 200170) 

NL number of points to be computed on each line 

NTIM when multiplied by time increment DELT represents the value of time on initial 
line 

NTOT number of time intervals to be computed 

KTPT if 1, results will not be punched on cards; if 2, results will be punched on cards 
for possible use in plotting 

Second input card (corresponding to read statement 200180) 


GAM 

ratio of specific heats 

GG 

acceleration due to gravity 

RR 

gas constant 

AJJ 

mechanical equivalent of heat 

PRS1 

reference pressure 

RHOl 

reference density 

AL0 

reference length 


Third input card (corresponding to read statement 200190) 

AKDS amplitude for — oscillation (cf. card 700220) 

Dr 
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I 

AKMK amplitude for Mach number oscillation at right boundary (cf. card 700210) 

AKPS amplitude for pressure oscillation at left boundary (cf . card 500140) 

AKS amplitude for entropy oscillation at left boundary (cf. card 500170) 

DSDTI initial value of 

Dr 

Fourth input card (corresponding to read statement 200200) 

FFDS frequency for — oscillation (cf. card 200350) 

Dr 

FFMK frequency for Mach number oscillation at right boundary (cf . card 200840) 

FFPR frequency for pressure oscillation at left boundary (cf. card 200340) 

FFS frequency for entropy oscillation at left boundary (cf. card 200360) 

DELT nondimensional time increment 

DELL distance between any two points on any base line; if properties on initial line 
are not constant, DELL = 0. 0 

EMI shock Mach number; if diffuser problem is not being run, EMI = 0. 0 
A21 diffuser exit to throat area ratio 

After the fourth input card, there are two options which are selected by setting DELL 
equal to either zero or to its value, the distance between two points on a base line. If 
DELL is zero, parameters at each location on the initial line are different. One card is 
required for each point. Hence, if the number of points on a base line (NL) is five, for 
example, five additional cards are required. 

So, for the first option: 

Input cards (corresponding to card 200230) 

Z(L) location of L^ 1 point on initial line 

it. 

R(L, 1) nondimensional acoustic velocity at L n point on initial line 

1U 

U(L, 1) nondimensional velocity at L U1 point on initial line 

ii. 

S(L, 1) nondimensional entropy at L n point on initial line 
TIM(l) nondimensional time on initial line 
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and for the second option: 


Fifth input card (corresponding to card 200260) 

AZ location of left boundary point, usually zero 

AR constant value of nondimensional acoustic velocity along initial line 

AU constant value of nondimensional velocity along initial line 

AS constant value of nondimensional entropy along initial line 

TIM(l) nondimensional time on initial line 


Fortran Program Listing 


c main program 

COMMON /CONST/ AO* AJJ, AKDS, AKPS, AKS* AMOK, CONI, C0N2, C0N3, 

1 CXI, DSDT , DSDTI, GAM, GG, GMM1 , GMP1, OMGD , OMGP, OMGS, 

2 PRS1 , PRSIL, RR, RH01, S3I, TIMO 

COMMON /PROP/ NIT (20,10) , R(20,10>, S(20,10), TIM(10), TIMR(IO), 
1 U ( 20 » 1 0 ) , Z ( 20 ) 

COMMON /RUN/ DELL, DELT, FFDS, FFPR , FFS , KTPT , NL , NT I M , NTOT 
COMMON /SUPSON/ AAST(lO), AKMK , AMKI , AXl(lO), A21, 

X Dl, DLS, DPP (10) , 

1 EMI, EM2 , EMXtlO), EMY(IO), FFMK , OMGM , RT , SPR 
DIMENSION PRSI20), RH0I20), TEMP ( 20 ) , ENT(IO), DST(IO), VEL(IO) 
DIMENSION AMK(IO), PLTD(840> 

1 CALL redy 
KDSTP= 1 
IBEG=1 

CALL BNET(NL,1) 

AAST ( 1 ) = AMK I 

21 IL5T= I BEG+? 

DO 24 I =2 , ILST 

T I M ( I >=NT0T 

TIMI I >=DELT*TIM{ I ) 

T I MR ( I )=TIM0*TIM( I ) 

NTOT =NT0T + 1 

IFf NT0T-NTIM)23,23,22 

22 I LST =1-1 
KDSTP = 2 
GO TO 28 

23 CONTINUE 

CALL BNET t NL , I ) 

DST ( I ) =DSDT 
AAST ( I )=AM0K 

24 CONTINUE 

28 I F ( ILST-1 129,87,29 

C PRINTED AND PUNCHED OUTPUT 

29 WR I TE ( 6 » 30 ) ( T I M ( I ) , I = I BEG, I LST ) 

30 F0RMAT(10H1 TIM = 8G15.6) 

WRITE(6*330) (TIMR(I) , I = I BEG , I LST ) 

330 FORMAT ( 10HJREL TIM= 8G15.6) 

WRITE (6, 340) ( DST ( I ) , I = I BEG , I LST ) 

340 FORMAT ( 10HJ DSDT = 8G15.6) 

I F ( EMI ) 32 , 36 , 32 


100020 

100030 

100040 

100050 

100060 

100070 

100080 

100090 

100100 

100110 

100120 

100130 

100140 

100150 

100160 

100170 

100180 

100190 

100200 

100210 

100220 

100230 

100240 

100250 

100260 

100270 

100280 

100290 

100300 

100310 

100320 

100330 

100340 

100350 

100360 

100370 

100380 

100390 

100400 

100410 

100420 
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in ii ■■ ■■ in 


in 


I llll I I III 



32 DO 34 I * I BEG » I LST 

34 AAST(I )=AST(AAST< I ) ,GMM1,GMP1 ) 

WRITE (6*350) ( AAST ( I ) , I = I BEG , I LST ) 

350 FORMAT ( 10HJ A/AST = 8G15.6) 

WRITE (6. 360) < AXl ( I ) , I = I BEG, I LST ) 

360 FORMAT ( 10HJ AX/A1 = 8G15.6) 

WRITE (6. 370) < DPP ( I ) . I-IBEG* I LST ) 

370 FORMAT < 10HJ DP/P * 8G15.6) 

WRITE (6. 380) ( EMX ( I ) . I *» I BEG. I LST ) 

380 FORMAT ( 10HJ MX = 8G15.6) 

WRITE (6. 390) ( EMY ( I ) ♦ I = I BEG . I LST ) 

390 FORMAT ( 10HJ MU = 8G15.6) 

36 GO TO (56.38) .KTPT 

38 K = 0 

DO 39 IMBEG.ILST 
K*K+1 

39 PLTD( K ) *T I M ( I ) 

I F ( EMI ) 41 » 56 *41 
41 DO 43 L=1.NL 

DO 43 I=IBEG.ILST 
K = K+1 

43 PLTD(K)=S(L.I) 

DO 45 I = I BEG . I LST 
K*K + 1 

PLTD ( K ) =AAST ( I ) 

K = K+1 

PLTD ( K ) = AX 1 ( I ) 

OK+1 

PLTD ( K ) =DPP ( I ) 

K = K+ 1 

45 PLTD ( K ) = EMX ( I ) 

56 DO 73 L=1 *NL 

WRITE (6,40) Z(L). ( R ( L . I ) , I = I BEG . I LST ) 

40 FORMAT! 1HKF7 . 4 , 2HR=8G1 5 . 6 ) 

WRITE (6,50) (U(L,I ) ,I=IBEG*ILST) 

50 FORMAT ( 10HJ U=8G15.6) 

WRITE (6,60) ( S ( L , I ) , I = IBEG, ILST) 

60 FORMAT ( 10HJ S=8G15.6) 

WR I TE ( 6 . 75 ) ( N I T ( L , I ) ,I = IBEG, ILST) 

75 FORMAT ( 10HJ NIT=8G15.6> 

DO 63 I = I BEG , I LST 
CN1=EXP(-GAM*S(L. I ) ) 

RHO ( I )=CN1*RH01* (R (L , I ) **C0N1 ) 

PRS ( I ) =CN1*PRS1*(R(L»I ) **C0N2 ) 

AA = AO*R (L , I ) 

TEMP ( I )=AA*AA/(GAM*GG*RR) 

ENT ( I ) =S ( L « I )*GAM*RR/AJJ 
VEL ( I ) =AO*U ( L , I ) 

AMK ( I ) =U ( L , I )/R(L,I ) 

63 CONTINUE 

WRITE16.65) ( ENT ( I) , I = I BEG , I L ST ) 

65 FORMAT ( 10HJ ENT=8G15,6) 

WRITE(6.70) (RHO( I)»I=IBEG»ILST1 
70 FORMAT ( 10HJ RH0=8G15,6> 

WR I TE ( 6 . 8 0 ) (PRS( I ) » I = I BEG » I LST ) 

80 FORMAT ( 10HJ PRS=8G15,6) 

WRITE (6. 90) ( TEMP ( I ) , I=IBEG, ILST ) 

90 FORMAT ( 10HJ TEMP=8G15.6) 

WRITE (6» 100) (VEL ( I ) , I = IBEG, ILST ) 

100 FORMAT ( 10HJ VEL=8G15.6) 

WRITE (6. 110) ( AMK ( I ) , I = I BEG, I LST ) 


100430 

100440 

100450 

100460 

100470 

100480 

100490 

100500 

100510 

100520 

100530 

100540 

100550 

100560 

100570 

100580 

100590 

100600 

100610 

100620 

100630 

100640 

100650 

100660 

100670 

100680 

100690 

100700 

100710 

100720 

100730 

100740 

100750 

100760 

100770 

100780 

100790 

100800 

100810 

100820 

100830 

100840 

100850 

100860 

100870 

100880 

100890 

100900 

100910 

100920 

100930 

100940 

100950 

100960 

100970 

100980 

100990 

101000 

101010 

101020 

101030 
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110 

F ORMAT ( 10HJ M0K = 8G15.6) 

101040 


R(L»1)=R(L,ILST) 

101050 


U(L,1)=U(L,ILST) 

101060 


S(L,1)=S(L,ILST) 

101070 


GO TO (73, 67), KTPT 

101080 

67 

DO 69 I = I BEG , I LST 

101090 


K = K+1 

101100 


PLTD ( K ) =TEMP ( I ) 

101110 


K = K+1 

101120 


PLTD ( K ) =RHO ( I ) 

101130 


K = K+1 

101140 


PLTD ( K ) =PRS.t I ) 

101150 


K = K+1 

101160 

69 

PLtD(K)=AMK( I ) 

101170 

73 

CONTINUE 

101180 


T IM ( 1 ) =T IM ( I LST ) 

101190 


GO TO (85,78) »KTP T 

101200 

78 

I F ( IBEG-1 >83,81,83 

101210 

81 

WR I TE ( 6 , 1 30 ) NL, K 

101220 

130 

FORMAT ( 2H$ 215) 

101230 

83 

CALL BCDUMP (PLTDt 1) ,PLTD‘(K) ,1 ) 

101240 

85 

I BEG=2 

101250 

87 

GO TO (21,88) ,KDSTP 

101260 

88 

GO TD (93,91) ,KTPT 

101270 

91 

WR I TE ( 6 , 1 30 ) K 

101280 

93 

GO TO 1 

101290 


END 

101300 


200020 

200030 

200040 

200050 

200060 

200070 

200080 

200090 

200100 

200110 

200120 

200130 

200140 

200150 

200160 

200170 

200180 

200190 

200200 

200210 

200220 

200230 

200240 

200250 

200260 

200270 

200280 

200290 

200300 

200310 

200320 


CONI. C0N2. CONS * 
GMP1 . OMGD* DMGP . OMGS. 

S t 20 * 1 0 ) » TIMUO). T I MR ( 1 0 ) * 

KTPT, NL, NT I M . NTOT 

A21. 

SPR 


1 


SUBROUTINE REDY 

COMMON /CONST/ AO. AJJ. AKDS . AKPS, AKS. AMOK. 

1 CXI. DSDT * DSDTI. GAM, GG, GMM1 , 

2 PRS1 » PRSIL, RR, RH01, S3 1 , TIMO 
COMMON /PROP/ NITI20.10). R(20»10>» 

1 U ( 20. 10 ) , Z ( 20 ) 

COMMON /RUN/ DELL, DELT. FFDS » FFPR » FFS 
COMMON /ALNE/ALIM.CA2 ,CA3 ,DIST ,SLP »ZLIM, 

1 CU( 19 ) ,CR ( 19) »CS( 19) »CURP ( 19 ) .CURMt 19) , 

2 ZL(20),RL(20),UL(20)»SL(20)»TL(20) 

COMMON /SUPSDN/ AAST(IO), AKMK , AMK I , AXltlO). 

X Dl, DLS, DPP ( 10), 

1 EMI, EM2, EMX(IC), EMY(IO), FFMK » OMGM » RT, 

PI = 3. 1415927 

INITIAL INPUT AND OUTPUT 
READ(5»10) NL.NTIM.NTOT.KTPT 
READ (5.20) GAM, GG, RR, 

READ t 5 , 20 ) AKDS, AKMK, . _____ 

READ! 5,20) FFDS, FFMK, FFPR, FFS, DELT, DELL, EMI, A21 
IF(DELL)5»3»5 
DO 4 L= 1 , NL 

READ (5,20) Z(L), R ( L , 

DELL=Z ( 2 ) -Z ( 1 ) 

GO TO 8 

READ(5»20) AZ, AR» AU, AS, TIM(l) 

CALCULATE INITIAL CONSTANTS 
DO 6 L= 1 ♦ NL 
CN1 =L-1 

Z ( L ) =AZ+CN1*DELL 
R (L ,1 )=AR 
U(L,1)=AU 


AJJ , 

PRS 1 

. RHDlt AL 

AKPS. 

AKS. 

DSDTI 

FFPR , 

FFS. 

DELT, DELL, 

) ♦ U(L 

.1 ) . 

S ( L , 1 ) ♦ TIM 
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c 

c 

c 


6 S ( L » 1 ) = AS 
8 DMGP=2»*PI*FFPR 
0MGD=2.*PI*FFDS 
0MGS=2.*PI*FFS 
R3 I L=R (1*1) 

S3 I =S ( 1.1 ) 

GMM1=GAM-1. 

GMP1=GAM+1. 

C0N1=2 • / GMM1 

C0N2=GAM#C0N1 

RHD0=RH01/GG 

A0=SQRT<GAM*PRS1*14A./RH00> 


TIMO=ALO/AO 
AMK I =U ( NL » 1 ) /R ( NL » 1 ) 

A AST (1)=AST(AMKI *GMM1 ,GMP1 ) 

WR I TE ( 6 » 1 00 ) NL * NT I M * NTOT * KTPT 
100 F ORMAT ( 8H1 NL = I3.5X.7HNTIM = I3*5X»7HNT0T 
WRI TE ( 6 * 1 10 ) DELL. DELT. PRS1, RH01, ALO 
110 FORMAT (8HKDELL = G16 . 8 , 3 X .7HDELT = G16 . 8 . 3X * 7HPRS 1 = G16.8.3X. 
17HRH01 = G16.8.3X.7H LO = G16.8) 

WRITEI6.120) GAM. GG » RR » AJJ 
120 F ORMAT ( 8HK GAM = G16.8.3X.7H G 
17H J = G16.8) 

WR I TE ( 6 . 130 ) DSDTI. AMK I . R3IL. S3I 
130 F0RMAT(8HKDSDTI= G1 6 . 8 » 3X . 7HAMK I = G16.8.3X.7HR3IL = G16, 


I3.5X.7HKTPT = 13) 


G16.8.3X.7H RR = G16.8.3X. 


17H S3 1 = G16.8) 

WR I TE ( 6 ♦ 140 ) AKDS » AKMK . AKPS. AKS 
1 AO FORMAT ( 8HKAKDS = G 16 . 8 ♦ 3 X ,7H AKMK = G1 6 . 8 , 3X . 7 HAKPS = G16.8.3X. 


17H AKS ' G16.8) 

WRITEI6.150) FFDS » FFMK » FFPR . FFS 
150 FORMAT ( 8HKFFDS = G 1 6 . 8 . 3X . 7HFFMK = G1 6 . 8 . 3X » 7HFFPR = G16.8.3X. 
17H FFS = G16.8) 

WRITE(6»160) AAST(l), AO, EMI. A21 
160 FDRMAT (8HKA/A* = G16.8.3X.7H AO = G16.8,/. 

1 8HK Ml = G16.8 ,3X,7HA2/A1= G16.8) 


CN1=EXP(-GAM*S3I ) 

PRSIL=CN1*PRS1*(R3IL**C0N2) 

CN1 = • 5*GMM 1 
C0N3=EXP ( CN1 ) 

CX 1 =CN 1 /GAM 

DIST = DISTANCE BETWEEN 2 POINTS IN THE GAS CIRCUIT 
DIST=DELL 

SLP = THE SLOPE OF THE BASE LINE IN THE GAS CIRCUIT 
SLP = 0 » 

CA2 = SORT ( 1 »+SLP**2 ) 

CA2 = 1 • 

CA3=CA2/DIST 
AL I M = Z ( 1 ) 

ZLIM=Z(NL) 

IF( EMI )21 ,25,21 

21 CALL SSBZ(GAM»GMM1,R(1,1 ) ,U( 1 .1 ) .5(1,1 ) ) 
0MGM=2.*PI*FFMK 
DSDT = DSDT I 
25 CONTINUE 
10 FORMAT (1615) 

20 FORMAT ( 8 E 1 0 » 0 ) 

RETURN 

END 


200330 
2 003 AO 
200350 
200360 
200370 
200380 
200390 
200A00 
200A10 
200A20 
200A30 
200AA0 
200A50 
200A60 
200470 
200A80 
200490 
200500 
200510 
200520 
200530 
2005A0 
200550 
200560 
200570 
200580 
200590 
200600 
200610 
200620 
200630 
2006A0 
200650 
200660 
200670 
200680 
200690 
200700 
200710 
200720 
200730 
2007A0 
200750 
200760 
200770 
200780 
200790 
200800 
200810 
200820 
200830 
2 008 AO 
200850 
200860 
200870 
200880 
200890 
200900 
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c 


c 


c 


SUBROUTINE BNET < NL , I > 

ORGANIZATION OF CALCULATIONS ON LINE CORRESPONDING TO TIME I 
COMMON /PROP/ NITI20.10). R<20»10), S(20*10). TIM(IO), TIMR(IO), 
1 U(20»10)» Z ( 2 0 ) 

COMMON/ALNE/ALIM, CA2 *CA3*DIST *SLP »ZL I M » 

1 CU( 19) ,CR< 19) *CS( 19) *CURP ( 19 ) .CURMt 19) * 

2 ZL(20) *RL(20) *UL(20) » SL ( 2 0 ) » TL ( 20 ) 

LF=NL-1 

I F ( 1-1)10,10.1 

1 J=I-1 

SET UP NET FOR LOCATIONS 1 
CALL GASPZ ( Z ( NL ) » TIM(I), 

1 NITtNL.I ) » T I MR ( I ) , NL, 

DO 8 M=2 ,LF 
L=LF+2-M 

CALL GASP (Z(L), TIMtl), R(L,I), 

8 CONTINUE 

CALL GASPA (Z(l), TIM(I), R t 1 , I ) , 

1 NIT< 1 , I ) , TIMR ( I ) , 1,1) 

SET UP BASE LINE CONSTANTS TO BE USED 
10 DO 28 LL=1,LF 
LR=LL+1 


TO NL ( F INAL LOCATION) 

R( NL, I ) , U ( NL , I ) , S ( NL » I ) 
I ) 


U ( L , I ) » S ( L , I ) , NITtL.I), L) 


U( 1 »I ) » S(1,I), 


IN CALCULATING THE NEXT LINE 


ZL(LL)=Z(LL) 

RL(LL)=R(LL,I ) 

UL (LL)=U(LL,I ) 

SL t LL ) =S ( LL » I ) 

TL ( LL) =TIM( I ) 

CU(LL)=U(LR,I ) -UL ( LL ) 

1 A CU(LL)=CA3*CU(LL) 

16 CR(LL)=R(LR,I )-RL(LL) 

22 CR(LL)=CA3*CR(LL) 

25 CURP(LL)=CU(LL)+CR(LL) 
CURMtLL)=CU(LL) -CR ( LL ) 
CS(LL)=CA3*(S(LR,I )-SL(LL) ) 
28 CONTINUE 

ZL ( NL ) =Z ( NL ) 

RL(NL)=R(NL,I ) 

UL ( NL ) =U ( NL , I ) 

SL t NL ) =S ( NL » I ) 

TL t NL ) =T I M ( I ) 

RETURN 

END 


300020 

300030 

300040 

300050 

300060 

300070 

300080 

300090 

300100 

300110 

300120 

300130 

300140 

300150 

300160 

300170 

300180 

300190 

300200 

300210 

300220 

300230 

300240 

300250 

300260 

300270 

300280 

300290 

300300 

300310 

300320 

300330 

300340 

300350 

300360 

300370 

300380 

300390 

300400 

300410 

300420 

300430 
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ti 

I 

l 


SUBROUTINE GASP < Z 3 , T 3 , R3 , U3 » S3 , I T , L ) 

C BASIC NET POINT CALCULATIONS 

COMMON /CONST/ AO. AJJ, AKDS . AKPS. AKS. AMOK. CONI, C0N2, C0N3, 

1 CXI. DSDT » DSDTI, GAM, GG, GMM1 , GMP1, OMGD, OMGP , OMGS. 

2 PRS 1 » PRSIL. RR » RHOl, S3I, TIMO 
COMMON /ALNE/AL I M, CA2 ,CA3 ,DIST ,SLP »ZLIM, 

1 CU( 19) ,CR t 19) »CSt 19) ,CURP ( 19 ) ,CURM( 19) , 

2 ZL( 20) ,RL ( 20) »UL( 20 ) ,SL(20 ) »TL( 20) 

CRIT=.001 

C INITIAL CALCULATIONS 

KDLN=1 
I T*0 
KDSTP= 1 
KLOP-1 
UP3=1000. 

RP3* 1000. 

TST2U=l»E+20 
TST2R“l«E+20 
U3=UL ( L ) 

R3“RL ( L ) 

C BEGINNING OF LOOP 

5 IT= IT+1 

C CALCULATIONS FOR POINT 1 

7 LB«L 
3 LA = LB- 1 

CN1«T3-TL(LA)+SLP*ZL(LA) 

CN2=U3+R3+UL ( LA ) +RL ( LA ) -CURP ( LA ) *ZL < LA ) 

AA»SLP*CURP(LA) 

8B=SLP*CN2-CNl*CURP(LA)-2. 

CC=2.*Z3-CN1*CN2 
Z 1 = QUAD (AA.BB.CC.-1..Z3 ) 

CN1=Z1-ZL(LA) 

T1=TL(LA)+SLP*CN1 
I F ( CN 1)12,19,19 
12 IF<ALIM-Z1)14,14,20 
14 LB=LB- 1 
GO TO 8 

19 U1=UL( LA)+CN1*CU( LA) 

R1=RL(LA)+CN1*CR(LA) 

S1=SL(LA)+CN1*CS(LA) 

GO TO 22 

20 T1 = T3-(Z3-ALIM)*(T3-U ) / ( Z3-Z 1 ) 

Z 1 = AL I M 

RELTM=T IM0*T1 

CALL GASPA(Z1,T1 ,R1 ,U1 ,S1 ,NDM , REL TM , LA ) 

KDLN=2 

C CALCULATIONS FOR POINT 2 

22 LB=L 
25 LC=LB+1 

CN1 =T3-TL ( LB)+SLP*ZL (LB ) 

CN2=U3-R3+UL(LB) _ RL(LB)-CURM(LB)#ZL(LB) 

AA=SLP*CURM( LB) 

BB=SLP*CN2-CNl*CURM(LB)-2. 

CC=2.*Z3-CN1*CN2 
Z2=QUAD( AA.BB.CC.l. ,Z3 ) 

CN 1 =Z 2-ZL < LB ) 

T2=TL(LB)+SLP*CN1 
I F ( Z2-ZL ( LC ) 132,32,27 

27 IF(Z2-ZLIM)2R,2R,34 

28 LB=LB+1 
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GO TO 25 

32 U2 = UL(I_B)+CN1*CU(LB) 

R2=RL(LB)+CN1*CR(LB) 

S2=SL(LB)+CN1*CS(LB> 

GO TO 37 

34 T2=T3-(Z3-ZLIM)*(T3-T2)/(Z3-Z2> 

Z2 = ZL I M 

CALL GASPZ (Z2»T2*R2*U2»S2*NDM*LC) 

KDLN = 2 

C CALCULATIONS FOR POINT 4 

37 GO TO (39,49) *KDLN 
39 LB=L 
41 LA=LB-1 

CN1 =U3+UL t LA ) -CU ( LA ) *ZL ( LA ) 

CN2=T3— TL(LA)+SLP*ZL(LA) 

AA=CU(LA)*SLP 
BB=CN1*SLP-CU(LA >*CN2-2. 

CC=2.*Z3-CN1*CN2 
Z4=QUAD( AA,BB*CC,-1. ,Z3) 

CN1=Z4-ZL ( LA ) 

I F ( CN1 ) 45 ,47,47 
45 LB=LB-1 
GO TO 41 

47 T4=TL(LA)+SLP*CN1 
U4=UL(LA)+CU(LA)*CN1 
S4=SL ( LA ) +CS ( LA ) *CN 1 
GO TO 51 

C ALTERNATE CALCULATIONS FOR POINT 4 WHEN BASE SLOPE DOES NOT = 0. 

49 D12=SQRT ( (T1-T2)**2+(Z1-Z2)**2) 

SL12=(T2-T1)/(Z2-Z1) 

CN1=SQRT ( 1.+SL12*SL12 )*(U2-U1 ) /D12 
CN2=T3-T2+SL12*Z2 
CN3=U3+U2-CN1*Z2 
AA=CN 1*SL 1 2 

BB=SL12*CN3-CNl*CN2-2. 

CC=2.*Z3-CN2*CN3 

Z4=QUAD(AA,BB,CC,-1.,Z3) 

CN2=Z2-Z4 

T4=T2-SL12*CN2 

U4=U2-CN1*CN2 

S4=S2-CN2*SQRT(1.+SL12*SL12)*(S2-S1)/D12 

KDLN=1 

C CALCULATE S3 

51 S3=S4+DSDT*(T3-T4) 

GO TO (54,57) » KLOP 
C CALCULATE U3 

54 P1=R1*2./GMM1+U1 
Q2=R2*2./GMM1-U2 
R13=.5* (R1+R3 ) 

P3 = Pl + Rl3#( GMM1 *DSDT * ( T3-T1 ) + (S3-Sl ) ) 

R2 3 = • 5* ( R2+R3 ) 

Q3 = Q2+R2 3*(GMM1*DSDT*( T3-T2 ) + (S3-S2) ) 

U3= • 5* ( P 3-Q3 ) 

KL0P=2 
GO TO 7 

C CALCULATE R3 

57 02=R2*2./GMM1-U2 
R2 3 = . 5* ( R2+R3 ) 

Q3=02+R23*(GMM1*DSDT*( T3-T2 )+ (S3-S2 ) ) 

KL0P=1 

R3=.5* (Q3+U3 ) *GMM 1 
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n r. 


I 



TEST FDR CONVERGENCE OF 

U3 AND R 3 

401230 


TST1U=ABS(U3-UP3 ) 


401240 


TST1R=ABS( R3— RP3 ) 


401250 


IF(TST1U-CRIT>62»65,65 


401260 

62 

I F ( TST 1R-CR I T)95»65»65 


401270 

65 

IF(TST2U-TST1U ) 7 1.71, 67 


401280 

67 

IF(TST2R-TST1R>73»73.84 


401290 

71 

WR I TE (6.600) Z3» T3 


401300 

600 

FORMAT ( 28HK* U IS NOT CONVERGING AT Z =G16 . 8 . 5 X » 2HT=G16 . 8 ) 

401310 


GO TO 75 


401320 

73 

WR I TE ( 6 »610 ) Z3, T3 


401330 

610 

FORMAT ( 28HK* R IS NOT CONVERGING AT Z=G16 . 8 » 5 X, 2HT=G16 .8 ) 

401340 

75 

CONTINUE 


401350 

76 

WR I TE ( 6 » 50 ) Zl, Rl, Ul, 

SI, T 1 

401360 


WRI TE ( 6 ,60 ) Z2, R2, U2, 

S2 , T 2 

401370 


WR I TE ( 6 .60 ) IT, Z4, U4 , 

S4, T 4 

401380 


WR I TE ( 6 »60 ) PI, 02. P3. 

Q3 

401390 


WRITEI6.60) Z3 » R3, U3» 

S3, T 3 

401400 

50 

FORMAT ( 4HL#* 5G20.8) 


401410 

60 

FORMAT ( 4H J 5G20.8) 


401420 

82 

GO TO ( 84,96 ) .KDSTP 


401430 

84 

UP3=U3 


401440 


RP3=R3 


401450 


T ST2U= TST 1U 


401460 


TST2R°TST1R 


401470 


GO TO 5 


401480 

95 

continue 


401490 

96 

IF1U3-R3) 101.101.98 


401500 

98 

WRITE(6»70)U3,R3,Z3,T3 


401510 

70 

FORMAT ( 30HL** STDP SUPERSONIC FLOW U=G 1 6 . 6 , 5X , 2HR=G 1 6 . 6 , 

401520 


15X.6HAT Z3=G16.6»5X»7HAND T3=G16.6) 

401530 


STOP 


401540 

101 

RETURN 


401550 


END 


401560 


SUBROUTINE GASPA ( Z 3 * T3 ,'R 3 , U3 . S3 * I T » REL TM ♦ L » I ) 

C LEFT BOUNDARY CALCULATIONS 

COMMON /CONST/ AO. AJJ. AKDS . AKPS* AKS. AMOK. CONI. C0N2. C0N3. 

1 CXI, DSDT , DSDTI, GAM, GG , GMM1, GMP1, OMGD, OMGP, OMGS, 

2 PR S 1 ♦ PRSIL, RR, RHOl, S3 I , TIMO 
COMMON/ALNE/ Al. IM.CA2.CA3.DIST ,SLP »ZLIM, 

1 CU( 19 ) ,CR ( 19) »CS t 19 ) »CURP ( 19 ) ,CURM ( 19 ) » 

2 ZL( 20) ,RL ( 20) ,UL ( 20 ) ,SL (20 ) »TL ( 20) 

CR I T= • 001 

INITIAL CALCULATIONS 
CALCULATIONS FOR S3 AND R3 
CN1=SIN(0MGP*RELTM) 

PRS3=PRSIL+CN1«AKPS*PRSIL 
CN2= ( PRS3/PRS1 )**CX1 
CN1=SIN(0MGS*RELTM) 

S3=S3I+AKS*CN1*S3I 
R3=CN2*(C0N3**S3) 

I F ( UL ( 1 ) -RL ( 1 ) ) A, 2,4 
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2 Z2=ZL ( 1 ) 500200 

T2 = TL ( 1 ) 500210 

U2=UL ( 1 ) 500220 

R2=RL ( 1 ) 500230 

S2=SL ( 1 ) 500240 

02=R2*2./GMM1-U2 500250 

R23 = . 5* ( R2+R3 ) 500260 

Q3 = Q2 + R23*tS3-S2 ) 500270 

U3 = 2 • *R3 /GMM 1-Q3 500280 

GO TO 25 500290 

4 I T=0 500300 

UP3=1000. 500310 

K D S T P = 1 500320 

TST2=l.E+20 500330 

U3=UL ( L ) 500340 

C BEGINNING OF LOOP 500350 

5 IT= IT+1 500360 

C CALCULATIONS FOR POINT 2 500370 

LB=L 500380 

8 LC=LB+1 500390 

CN1=T3-TL(LB)+SLP*ZL(LB) 500400 

CN2=U3-R3 + UL ( LB) -RL ( LB l-CURM ( LB ) *ZL ( LB 5 00410 

AA=SLP*CURM( LB) 500420 

BB = SLP*CN2-CN1*CURM< LB) -2. 500430 

CC=2.*Z3-CN1*CN2 500440 

Z2=QUAD(AA,BB»CC»1.»Z3) 500450 

IF(Z2-ZL(LC) ) 12*12.10 500460 

10 LB=LB+1 500470 

GO TO 8 500480 

12 CN1 =Z 2-ZL ( LB ) 500490 

T2=TL(LB)+SLP*CN1 500500 

U2=UL(LB)+CN1*CU(LB) 500510 

R2=RL( LB)+CN1*CR ( LB) 500520 

S2=SL(LB)+CN1*CS(LB) 500530 

Q2=R2*2./GMM1-U2 500540 

C CALCULATE U3 500550 

R23=.5*(R2+R3) 500560 

Q3=Q2+R23*(GMM1*DSDT*(T3-T2)+(S3-S2) ) 500570 

U3=2.*R3/GMM1-Q3 500580 

C TEST FOR CONVERGENCE OF U3 500590 

TST1=ABS(U3-UP3) 500600 

I F ( TST1-CRIT >25.15.15 500610 

15 I F ( TST2-TST1 ) 16.17,17 500620 

16 WR I TE ( 6 » 50 ) Z2, R2, U2» S2» T2 500630 

WR I TE ( 6 .60 ) IT, Q2, Q3 500640 

WR I TE ( 6 , 6 0 ) Z3 , R3, U3, S3, T3 500650 

50 F ORMAT ( 4HL** 5G20.8) 500660 

60 F ORMAT ( 4H J 5G20.8) 500670 

GO TO ( 17 . 26 ) .KDSTP 500680 

17 UP3 =U3 500690 

TST2=TST1 500700 

GO TO 5 500710 

25 IF(U3-R3)34,34,33 500720 

33 U3=R3 500730 

34 CONTINUE 500740 

26 RETURN 500750 

END 500760 
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SUBROUTINE' GASPA t Z 3 * T3 .R3 ,U3 , S3 » I T ,RELTM ,L . I ) 

C LEFT BOUNDARY SUBROUTINE FOR THE SUPERSONIC BUZZ PROGRAM 

COMMON /CONST/ AO» AJ J * AKDS * AKPS, AKS. AMOK. CONI. CON2. C0N3. 

1 CXI. DSDT . DSDTI, GAM, GG, GMM1 , GMP1, OMGD, OMGP , OMGS. 

2 PRS1 » PRSIL, RR, RHOl, S3 1 , TIMO 
C0MM0 N/AlNE/ALIM»CA2.CA3,DIST,SLP.ZLIM, 

1 CU< 19) ,CR< 19) .CS( 19) »CURP ( 19 ) »CURM( 19) , 

2 ZL ( 20 ) »RL(20) »UL(20)»SL(20) ,TL(20) 

COMMON /SUPSON/ AAST(IO), AKMK, AMK I , AXl(lO), A21, 

X Dl, DLS, DPP ( 10 ) , 

1 EMI, EM2 « EMX(IO). EMY(IO), FFMK , OMGM, RT, SPR 
CRIT=.00l 

C INITIAL CALCULATIONS 

CNN=GMPl/GMMl 
I F ( UL t 1 ) -RL ( 1 ) )4»2,4 
2 Z2»ZL<1) 

T2=TL ( 1 ) 

U2 = UL ( 1 > 

R2=RL ( 1 ) 

S2*SL ( 1 ) 

Q2=R2*2./GMM1-U2 

R23=.5*<R2+R3) 

Q3=Q2+R23*(S3-S2 ) 

U3=2.*R3/GMM1-Q3 
GO TO 25 
A I T = 0 
KDSTP= 1 
TST2= 1 • E+20 

Q3=2.*RL(L)/GMM1-UL(L) 

YY 1=0. 

XXI = .99*03 
QP3=Q3 

C BEGINNING OF LOOP 

5 I T= I T+l 
XX2 = 03 

C CALCULATIONS FOR R.U, AND S AT DIFFUSER EXIT 

R3 = Q3 + S0RT <CNN*RT*RT-.5*GMM1*Q3*Q3 ) 

R3=R3/CNN 

U3=2.*R3/GMM1-Q3 

EM3=U3/R3 

D3= 1 • + . 5*GMM1*EM3*EM3 
D3 = EM3/(D3**( .5*CNN) ) 

DSR=ALOG( A21*D3/D1 ) 

DLS=DSR/GAM 

S3=SPR+DLS 

C CALCULATIONS FOR POINT 2 

LB=L 

8 LC=LB+1 

CN1=T3-TL (LB)+SLP*ZL (LB ) 

CN2=U3-R3+UL ( LB) -RL ( LB)-CURM ( LB)*ZL (LB ) 

AA=SLP*CURM(LB) 

BB=SLP*CN2-CNl*CURM(LB)-2. 

CC=2.*Z3-CN1*CN2 

Z2 = QUAD( AA.BB.CC.l. , Z 3 ) 

I F ( Z2-ZL ( LC ) >12,12,10 
10 LB=LB+1 
GO TO 8 

12 CN1 =Z 2-ZL ( LB ) 

T2=TL(LB)+SLP*CN1 

U2=UL(LB)+CN1*CU(LB) 

R2=RL(LB)+CN1*CR(LB) 

S2=SL(LB)+CN1*CS(LB> 

Q2=R2*2, /GMM1-U2 
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R23= . 5* ( R2+R3 ) 

600650 



Q3 = Q2 + R23*(GMM1*DSDT*( T3-T2 )+ (S3-S2 ) ) 

600660 

c 


TEST ON CONVERGENCE OF THE Q COMPATIBILITY RELATION 

600670 



YY2=Q3-QP3 

600680 



TSTl = ABSt YY2/Q3 ) 

600690 



IF(TST1-CRIT)25,17,17 

600700 


17 

QP3=Q3 

600710 



Q3=(XX1*YY2-XX2*YY1 )/( YY2-YY 1 ) 

600720 



YY1=YY2 

600730 



XXI =XX2 

600740 



TST2=TST 1 

600750 



GO TO. 5 

600760 


25 

IF(U3-R3)34»3A,33 

600770 


33 

U3 = R3 

600780 


34 

CONTINUE 

600790 


26 

EMX ( I )=FMX(EM1 .DSR.GAM.GMM1.GMP1 ) 

600800 



CN1 =EMX ( I )*EMX( I ) 

600810 



DX=1.+.5»GMM1*CN1 

600820 



DX=EMX ( I ) / ( DX** ( • 5*CNN ) ) 

600830 

c 


CALCULATION FOR SHOCK LOCATION 

600840 



AX 1 ( I )=D1/DX 

600850 



I F t 1.-AX1 ( I ) ) A6 .46,42 

600860 


42 

WRITEI6.70) 

600870 


70 

FORMAT 1 26HL** AX/A1 IS LESS THAN 1.) 

600880 


44 

STOP 

600890 


46 

I F ( AX 1 ( I ) -A2 1 ) 49 » 49 . 47 

600900 


47 

WRITEI6.80) 

600910 


80 

FORMAT ( 32HL** AX/A1 IS GREATER THAN A2/A1 ) 

600920 



GO TO 44 

600930 


49 

EMXST=SQRT (GMPl#CNl/(2.+GMMl*CNl ) ) 

600940 



EMYST=1./EMXST 

600950 



CN1=EMYST*EMYST 

600960 



EMY t I ) =SQRT(2.*CNl/( GMP 1-GMM 1 *CN 1 ) ) 

600970 

c 


CALCULATION FOR PRESSURE RECOVERY 

600980 



DPP ( I ) =A21*D3/D1-1. 

600990 



RETURN 

601000 



END 

601010 


c 


c 


SUBROUTINE GASPZ < Z 3 , T3 » R 3 , U3 , S3 ♦ I T , REL TM , L , I ) 

RIGHT BOUNDARY CALCULATIONS 

COMMON /CONST/ AO» A J J ♦ AKDS ♦ AKPS* AKS, AMOK ♦ CONI. CON2. C0N3, 

1 CXI. DSDT . DSDTI, GAM, GG, GMM1 , GMP1 , OMGD, OMGP, OMGS. 

2 PRS 1 ♦ PRSIL, RR, RHOl, S3 I » TIMO 
C0MM0N/ALNE/ALIM.CA2.CA3 ,DIST .SLP .ZLIM, 

1 CU( 19) ,CR ( 19) ,CS( 19) .CURPI 19) .CURMt 19) » 

2 ZLI20) »RL(20) ♦ UL t 20 ) »SL(20) * T L ( 2 0 ) 

COMMON /SUPSON/ AAST(IO). AKMK , AMK I , AXl(lO). A21, 

X Dl, DLS, DPP (10), 

1 EMI, EM2, EMX(IO), EMY(IO), FFMK, OMGM, RT , SPR 
II U3 IS NOT KNOWN, FIND U3, R3, AND S3 
CRIT=»00l 

INITIAL CALCULATIONS 
I T= 0 

RP3=1000. 

KDSTP= 1 

LM=L-1 

TST2=l.E+20 

AMOK = AMKl+AKMK*SIN< OMGM*RELTM ) 

DSDT = DSDTI+AKDS*SIN< OMGD*RELTM) 
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R3=RL(L> 

700230 


U3=AM0K*R3 

700240 


BEGINNING OF LOOP 

700250 

5 

I T= I T+ 1 

700260 


CALCULATIONS FOR POINT 4 

700270 


LB = L 

700280 

7 

LA=LB-1 

700290 


CN1=U3+UL(LA)-CU<LA)*ZL<LA> 

700300 


CN2=T3-TL(LA)+SLP*ZL(LA) 

700310 


AA=CU(LA)*SLP 

700320 


BB=CN1*SLP-CU(LA )*CN2-2. 

700330 


CC=2.*Z3-CN1*CN2 

700340 


Z 4= QUAD (AA*BB»CC»~1 . *Z3 ) 

700350 


IF(ZL(LA)-Z4)9»9»8 

700360 

8 

LB=LB-1 

700370 


GO TO 7 

700380 

9 

CN1=Z4-ZL (LA) 

700390 


T4=TL(LA)+SLP*CN1 

700400 


U4=UL(LA)+CU( LA)*CN1 

700410 


S4=SL(LA)+CS(LA)*CN1 

700420 


S3=S4+DSDT# ( T3-T4 ) 

700430 


CALCULATIONS FOR POINT 1 

700440 


LB=L 

700450 

11 

LA=LB-1 

700460 


CN1=T3-TL ( LA ) +SLP*ZL ( LA ) 

700470 


CN2=U3+R3+UL(LA)+RL(LA)-CURP( LA)*ZL(LA) 

700480 


AA=SLP*CURP ( LA) 

700490 


BB=SLP*CN2-CN1*CURP ( LA ) -2 • 

700500 


CC=2.*Z3-CN1*CN2 

700510 


Z1=QUAD< AA,BB»CC»-1. ,Z3 ) 

700520 


CN1=Z1-ZL (LA) 

700530 


IFICNl >12*13*13 

700540 

12 

LB=LB-1 

700550 


GO TD 11 

700560 

13 

T1=TL(LA)+SLP*CN1 

700570 


Ul=UL ( LA ) +CN1*CU ( LA ) 

700580 


R 1 =RL ( LA) +CN1*CR ( LA ) 

700590 


S1=SL(LA)+CN1*CS(LA) 

700600 


Pl = Rl*2./GNiMl+ul 

700610 


CALCULATIONS FOR U3 AND R3 

700620 


R13=.5*(R1+R3) 

700630 


P3=P1+R13* (GMM1*DSDT* ( T3-T1 ) + ( S3-S 1 ) ) 

700640 


U3=2 • / (GMM1*AM0K) 

700650 


U3=P3/ ( 1 ,+U3 ) 

700660 


R3=U3/AM0K 

700670 


TEST FOR CONVERGENCE OF R3 

700680 


TST1=ABS(R3-RP3) 

700690 


IF(TST1-CRIT)25.15,15 

700700 

15 

IF( TST2-TST1 ) 16.17,17 

700710 

16 

WRITEI6.50) Zl, Rl. U1 , SI. T1 

700720 


WR I TE ( 6 » 6 0 ) IT, Z4, U4 , S4, T4 

700730 


WRI TE ( 6 .60 ) PI, P 3 

700740 


WR I T E ( 6 ♦ 6 0 ) Z 3 , R3, U3 , S3, T3 

700750 

50 

F ORMAT ( 4HL** 5G20.8) 

700760 

60 

F ORMAT ( 4H J 5G2C.8) 

700770 


GO TO ( 17,26) .KDSTP 

700780 

17 

RP3=R3 

700790 


TST2=TST1 

700800 


GO TO 5 

700810 

25 

CONTINUE 

700820 

26 

RETURN 

700830 


END 

700840 
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SUBROUTINE SSBZ ( GAM.GMM1 ,R2 »U2 »S2 ) 

INITIAL CALCULATIONS FOR THE SUPERSONIC DIFFUSER PROBLEM 
COMMON /SUPSON/ AAST(IO)* AKMK, AMKI , AXl(lO), A21, 

X Dl. DLS, DPP( 10) , 

1 EMI, EM2 » EMX(IO), EMY(IO), FFMK , OMGM, RT , SPR 
GMP1=GAM+1, 

EM2=U2/R2 

PHE=R2*R2*(2./GMM1+EM2*EM2) 

RT=SQRT ( • 5*GMM1*PHE ) 

CN1 *• 5*GMP1 /GMM1 

Dl = 1 •+• 5*GMM 1*EM1 *EM1 

D1=EM1/(D1**CN1 ) 

D2= 1 • + • 5*GMM 1*EM2 *EM2 
D2=EM2/<D2**CN1 ) 

DSR=ALOG( A21*D2/D1 ) 

DLS=DSR/GAM 

SPR=S2-DLS 


10 


20 


GMX=.5*(1.+EM1) 

EMX(1)=FMX(GMX,DSR ,GAM ,GMM1 , GMP1 ) 
CN2=EMX<1)*EMX(1) 

DX= 1 •+• 5*GMM1*CN2 
DX*EMX( 1 ) /<DX**CN1 ) 

AXl (1)=D1/DX 

EMXST=S0RT(GMP1*CN2/ ( 2.+GMMl*CN2 ) ) 

EMYST=1./EMXST 

CN2=EMYST*EMYST 

EMY ( 1 ) =SQRT ( 2.*CN2/ { GMP 1-GMM 1*CN2 ) ) 
DPP(1)=A21*D2/D1-1. 

WR I TE ( 6 * 10 ) 

FORMAT ( A8HL INITIAL OUTPUT FOR THE SUPERSONIC 


BUZZ PROGRAM) 


WR I TE ( 6 * 20 ) AXl ( 1 ) , EM2 , EMX ( 1 ) * 
F0RMAT(8HJAX/A1= G16.8*3 x *7H M2 


EMY ( 1 ) 

= G16.8,3X,7H 


MX = G16.8.3X, 


17H MY = G16.8) 
WR I TE ( 6 • 30 ) DLS, 
30 FORMAT ( 8HK DLS = 
17HDP/P = G16.8) 
RETURN 
END 


RT, SPR, DPP ( 1 ) 
G16.8,3X,7H RT 


= G16.8,3X,7H SPR = G16.8.3X, 
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FUNCTION FMX(GMX»DSR,GAM,GMM1 ,GMP1 ) 900020 

C CALCULATION OF SHOCK MACH NUMBER FROM A KNOWN ENTROPY RISE ( DSR ) 900030 

CRIT=.0001 900040 

I T = 0 900050 

EMXPR= 1 • E+20 900060 

EMX=GMX 900070 

CN1 =GAM/GMM1 900080 

CN2=GMMl/GMPl 900090 

1 I T= I T+l 900100 

CN3=EMX*EMX 900110 

F=AL0G(2./(GMP1*CN3)+CN2) 900120 

F=AL0G(2.*GAM*CN3/GMP1-CN2 ) /GMM1+CN1*F-DSR 900130 

FPR=EMX/ t 2.*GAM»CN3-GMM1 ) 900140 

FPR=FPR-l./(EMX*(2.+GMMl*CN3) ) 900150 

FPR=4.*CN1*FPR 900160 

TEST=ABS( EMX-EMXPR) - 900170 

IF(TEST-CRIT)5,5*4 900180 

4 EMXPR=EMX 900190 

EMX=EMX~F/FPR 900200 

GO TO 1 900210 

5 FMX=EMX 900220 

RETURN 900230 

END 900240 


FUNCTION AST(AMK*GMM1,GMP1) 1000020 

C CALCULATION FOR CRITICAL AREA RATIO 1000030 

AST=.5*AMK*AMK*GMM1+1. 1000040 

AST=2.*AST/GMP1 1000050 

CN1 = .5*GMP1-/GMM1 1000060 

AST=(AST**CN1)/AMK 1000070 

RETURN 1000080 

END 1000090 



FUNCTION QUAD(AA,BB,CC,SGN,CLOS) 

1100020 


SOLUTION FOR QUADRATIC EQUATION 

1100030 


I F ( AA ) 8 » 5 * 8 

1100040 

5 

QUAD=-CC/B6 

1100050 


GO TO 20 

1100060 

8 

WRK=2.*AA 

1100070 


DISC=SQRT (BB*BB-4.*AA*CC ) /WRK 

1100080 


WRK=-BB/WRK 

1100090 


X 1 = WRK + D I SC 

1100100 


X2=WRK-DI SC 

1100110 


TST1=SGN*(X1-CL0S) 

1100120 


I F t TST 1 >16.10,10 

1100130 

10 

TST2 = SGN* ( X2-CL0S ) 

1100140 

. 

IFITST2) 14,12.12 

1100150 

12 

IFITST1-TST2) 14,14,16 

1100160 

14 

QUAD= X 1 

1100170 


GO TO 20 

1100180 

16 

QUAD=X2 

1100190 

20 

RETURN 

1100200 


END 

1100210 
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V u 


APPENDIX D 


COMPUTER PROGRAM FOR NONSTEADY LIQUID FLOW 

The following is the computer program used for the example for liquid dynamics. It 
is given for reference purposes only. This program may be compared to the program in 
reference 3. 


COMMON /CHIN/ AK* AKK (11)* CNKH * LBEG * LLST, LL2 » I NC CIO)* 

1 KDSTP. KTQH * NSAV. NTOT* TIM(IO). TMPD 
COMMON /RUNWH/ AKH(20>, AKAPT(20), AREA ( 10 ) * DELT ♦ 

1 HL ( 1 1 » 20 ) * HR (11*20) » KST(ll), KTRL ( 10 ) » N I NT » NPTS. 

2 00(11*20) * SX(10). SY ( 10 ! * TMEI20) 

EQUIVALENCE (KT1»KTRL(1 ) ) , ( KT2 *KTRL ( 2 ) ) * ( KT3 »KTRL ( 3 ) ) » 

1 (KT4.KTRLI4) ) , ( KT5 »KTRL ( 5 ) ) 

5 CALL WHIN 

CALCULATE CONTROLS TO DETERMINE WHEN INITIAL WAVE ARRIVES AT EACH 
LOCATION 
DO 14 1=1,5 
KST ( I ) = 1 

I F ( KT5 -I >11,8,10 
8 KST ( I ) =2 
GO TO 14 

10 KIN= I 
KFIN=KT5-1 
GO TO 12 

11 KIN=KT5 
KF I N= I -1 

12 DO 13 K=K I N , KF I N 

13 KST ( I ) =KST <11 + 1 NC < K ) 

14 CONTINUE 
KDSTP= 1 
NTO T = 0 
LBEG= 1 

18 LLST = LBEG+7 
CALL WHAM 

GO TO (38,38,48) , KDSTP 
38 CALL WHOUT ( KDSTP ) 

GO TO (45,48) , KDSTP 
45 LBEG=LL2 
GO TO 18 
48 GO TO 5 
END 
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SUBROUTINE WHIN 

COMMON / CNRD/ ALEN(IO), ANP, BA* BO* CONI* 

1 FF. GG» PI * QGG , TAUN 

COMMON /CHIN/ AK » AKK(11>* CNKH * LBEG , LLST * LL2 , INCflO), 
1 KDSTP * KTQH. NSAV. NTOT* TIM(IO)* TMPD 
COMMON /RUNWH/ AKHI20), AKAPT ( 20 ) . AREA ( 10 ) » DELT. 

1 HL ( 1 1 . 20 ) * HR ( 11*20) * KST(ll), KTRL<10), NINT. NPTS* 

2 QQ ( 11,20) ♦ SX(IO), SY(10), TMEI20) 

DIMENSION D I A ( 1 0 ) , TH(lO), VCW(IO), VW(lO) 

EQUIVALENCE (KT1 ,KTRL ( 1 ) > » ( KT2 ,KTRL < 2 ) > , ( KT3 » KTRL ( 3 ) ) * 

1 ( KT4 *KTRL ( A ) ) , ( KT 5 , KTRL < 5 ) ) 

EQUIVALENCE (HA,HTANK> 

C READ INPUT 

1000 READ (5,100) NINT, NPTS, KTRL 
100 FORMAT (1615) 

N I =N I NT+1 

READ (5*200) DELT, BO, BA, FF, CONI, TAUN, ANP 
READ ( 5 ,200 ) GG, QGG, HTANK, El, RO 
200 F0RMAT(8E10.0) 

READ! 5,200) AK, ( AKK (I) * I = 1 * N I ) 

READ(5,200) ( ALEN ( I ) , I = 1 ,N I NT ) 

PI=3. 1415927 
KTQH= 1 

I F ( ANP 11008,1005,1005 
1005 PERD=1 ,/FF 

TMPD=ANP*PERD 
GO TO 1009 

1008 TMPD=NPTS+1 
TMPD=TMPD*DELT 

1009 CNKH = AKK(NI ) 

GO TO (1,3) , KT3 

1 READ ( 5 , 200 ) ( D I A ( I ) , I = 1 , N I NT ) 

RFAD ( 5 , 200 ) ( TH ( I) , I = 1 , N I NT ) 

C INITIAL CALCULATIONS 

VO = SORT ( GG*AK*144»/R0) 

DO 2 1 = 1 , N I N T 

AREA! I )=PI *P I A ( I ) *D I A ( I )/576» 

VCW ( I ) =SQRT ( 1 ,+DI A( I ) *AK/ ( TH ( I ) *F1 I ) 

2 VW ( I ) =VO/VCW ( I ) 

GO TO 4 

3 READ(5,200) ( VW ( I ) , I = 1 , N I NT ) 

READ(5,200) ( ARE A ( I ) , I = 1 , N I N T ) 

4 NSAV = 0 

DO 7 1=1, NINT 

SX ( I ) = VW ( I ) / ( C-C*ARrA ( I ) ! 

SY ( I ) =-SX ( I ) 

T I M ( I )=ALTN< I ) /VW ( I i 
INC ( I ) =T IM ( I ) /DELT+.5 
IF ( NSAV-I NC ( 1 ) ) 6 , 7,7 

6 NSAV= INC ( I ) 

7 CONTINUE 
LL2=NSAV+2 

GO TO ( 17,20,20) ,KT4 
17 AKAor ( 1 ) =B0 

AKK ( 1 ) =C0N1* AKAf'T ( 1 i 
20 HR ( i , 1 ) =HA 
i .'*! = 1 

DO 23 1=1 , N I 
HL(I,i) = HR(IN’,l) 

I M- I 
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HR ( I»1)=HL(IM*1)-(QGG/AKK(I ) ) **2 
23 CONTINUE 
C INITIAL OUTPUT 

WR I TE ( 6 »210 ) NINT, NPTS, (KTRL( I ) *1 = 1*5) 

WR I TE ( 6 »220 1 DELT* BO, BA, FF, CONI, TAUN, ANP 
WR I TE ( 6 ,230 ) GG , QGG, HA, HRlNI.l), EX, RO , AK 
WRITE! 6*240 ) < AKK ( I ) , I = 1 ,NI ) 

210 FORMAT (8H1NINT = I3,5X,7HNPTS = 13,/, 

18HKKTRL = 13,917) 

220 F0RMAT(8HKDELT = G13.6.5X.7H BO = G13.6,5X,7H BA = G13.6,5X» 
17H FF *= G13.6,/, 

28HKC0N1 = G13.6,5X,7HTAUN = G13.6,5X,7H ANP = G13.6) 

230 FORMAT ( 8HK G - G13.6,5X*7H QGG = G13 .6 , 5X , 7HHTANK= G13.6.5X, 
17H HC = G13.6,/, 

28HK El = G13.6,5X,7H RO = G13.6,5X,7H K = G13.6) 

240 FORMAT ( 8HK KA = G13.6,5X,7H KB = G13.6,5X,7H KDE = G13.6,5X* 
17H KG = G13.6,5X,7H KH = G13.6) 

GO TO (43,48) ,KT3 
43 WR I TE < 6 ,235 ) VO 

WR I TE ( 6 *245 ) ( DI A ( I ) » I«1 ,NI NT ) 

WRI TE ( 6 *260 ) ( TH ( I ) , I = 1 »N I NT ) 

WRITE<6,280) «VCW( I ) ♦ I»1 ,NINT ) 

48 WRITE(6,250) < ALEN (I) , I =1 »N I NT ) 

WR I TE ( 6 » 270 ) ( AREA ( I ) , I = 1 ,N I NT ) 

WR I TE ( 6 *290 ) ( VW < I ) , I =1 »N I NT ) 

WR I TE ( 6 * 300 ) <SX< I ) ,1=1, NINT ) 

WRITE(6,310) <SY< I ) ,1 = 1, NINT) 

WR I TE ( 6 , 320 ) < T I M ( I ) , I = 1 ,N I NT ) 

WR I TE ( 6 , 330 ) ( INC ( I ) , 1 = 1 ,NINT 1 

235 FORMAT ( 8HK VO = G16.8,//, 

22HK ,14X,1H1,21X,1H2,21X,1H3,21X,1H4) 

245 FORMAT ( 8H J D = G16.8 ,4G22.8 ) 

260 FORMAT ( 8HK TH = G16. 8 ,4G22 . 8 ) 

280 FORMAT ( 8HK VCW = G16, 8 ,4G22 . 8 ) 

250 FORMAT ( 8HK L = G16. 8 ,4G22 .8 ) 

270 FORMAT ( 8HK A = G16.8 ,4G22 .8 ) 

290 FORMAT ( 8HK VW = G16 . 8 ,4G22 . 8 ) 

300 FORMAT ( 8HK SX = G16. 8 ,4G22 . 8 ) 

310 FORMAT ( 8HK SY = G16.8 ,4G22.8 ) 

320 FORMAT ( 8HK TIM = G16 . 8 ,4G22 . 8 ) 

330 FORMAT ( 8HK INC = G16.8 ,4G22 . 8 ) 

RETURN 

END 
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SUBROUTINE WHAM 300010 

C CHARACTERISTIC SOLUTION OF WATER HAMMER BY LATTICE POINT METHOD 300020 

COMMON /CNRD/ ALEN(IO), ANP, BA, BO, CONI, 300030 

1 FF, GG, PI, QGG, TAUN 300040 

COMMON /CHIN/ AK , AKK(ll), CNKH , LBEG , LLST , LL2 , I NC (10), 300050 

1 KDSTP, KTQH , NSAV , NTOT * TIM(IO), TMPD 300060 

COMMON /RUNWH/ AKHI20), AKAPT(20), AREA(IO), DELT, 300070 

1 HL ( 1 1 , 20 ) , HR ( 11,20) , KST(ll), KTRL(IO), NINT, NPTS, 300080 

2 QQ (11,20), SX(10), SY(10), TME(20) . 300090 

EQUIVALENCE (AKA, AKK ( 1 ) ) 300100 

EQUIVALENCE ( KT 1 , KTRL ( 1 ) ) , ( KT2 ,KTRL ( 2 ) ) , ( K T 3 , KTRL ( 3 ) ) , 300110 

1 ( KT4 *KTRL ( 4 ) ) , ( KT5 ,KTRL ( 5 ) ) 300120 

DIMENSION CX(ll), CY(ll) 300130 

NI=NINT+1 300140 
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n n 


c beginning of loop on time 

DO 90 lT=LBEG*LLST 
NT0T=NT0T+1 
IF(NT0T-NPTS)4,4»2 

2 LLST=LT-1 
KDSTP = 2 

IF(LLST-LBEG)3,95»95 

3 KDSTP=3 
GO TO 95 

4 AT=NT0T-1 
T I ME = AT *DELT 
TME.( LT ) X T I ME 
AKH ( LT ) = AKK ( N I ) 

AKAPT<LT)=AKA/C0N1 
I F t T I ME-TMPD ) 6.20,20 

6 GO TO ( 7, 8,20) ,KT4 

7 WRK=2.*PI*FF*TIME 
AKAPT ( LT ) =BO+BA*S IN ( WRK ) 

AKA=CONl*AKAPT(LT) 

GD TO 20 

8 AKH(LT)=CNKH*(1.-AT/TAUN) 

IF< AKH(LT) ) 9 , 9 , 1 2 

9 AKH ( LT ) =0, 

KTQH=2 

12 AKKINI ) = AKH t LT ) 

BEGINNING OF LOOP TO CALCULATE EACH LOCATION AT A GIVEN TIME 
(TME(LT) ) 

20 DO 85 1 = 1 , N I 
I F ( I -N I NT ) 21 , 2 1 , 2 3 

21 LPST=MAXO(LT-INC( I ) ,1 ) 

CY ( 1 )=HH 1+1 ,LPST )+SY( I >*QQ( 1+1 ,LPST ) 

CX( 1 + 1 ) =HR ( I ,LPST )+SX( I )*QQ< I ,LPST ) 

23 IF(KST< I l-NTOT >27,27,24 

24 QQt I »LT ) =QGG 
HL ( I »LT ) =HL ( I ,1 ) 

HR ( I ♦ LT ) =HR t 1,1) 

GO TO 85 

27 GO TO (32,42,52,62,72) , I 
32 TRM=AKK ( 1 ) *S Y ( 1 ) 

WRK=TRM*TRM-4.* (CY ( 1 ) -HL ( 1 ,1 ) ) 

QQ( 1 ,LT )= . 5*AKK ( 1 )* ( TRM+SORT ( WRK ) ) 

HR ( 1 ,LT ) =CY ( 1 ) -QQ< 1 , LT ) *SY ( 1 ) 

GO TO 85 

42 XX=CX ( 2 ) -CY ( 2 ) 

YY = SX ( 1 ) -SY ( 2 ) 

QQ ( 2 ,LT ) =OCALC ( AKK ( 2 ) , XX , YY ) 

HL(2»LT>=CX(2)-SX(1>*QQ(2,LT) 

HR(2,LT)=CY(2)-SY(2) *QQ ( 2 , LT ) 

GO TO 85 

52 XX=CX ( 3 ) -CY ( 3 ) 

YY=SX ( 2 ) -SY ( 3 ) 

QQ( 3»LT)=QCALC (AKK( 3) ,XX,YY) 

HL(3,LT)=CX(3)-SX(2) *QQ ( 3 ,LT ) 

HR( 3,LT)=CY(3)-SY (3) *QQ ( 3 ,LT ) 

GO TO 85 

62 XX=CX ( 4 ) -CY ( 4 ) 

Y Y = SX ( 3 > -SY ( 4 ) 

QQ (4,LT)=0CALC(AKK(4) ,XX,YY) 

HL(4,LT)=CX(4)-SX(3)*QQ(4,LT) 

HR(4,LT)=CY(4)-SY(4)*QQ(4,LT) 

GO TO 85 

72 GO TO ( 78,76) »KTQH 
76 QQ ( 5 , L T ) = 0 . 

GO TO 79 
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78 XX = CX ( 5 ) -HR (5*1) 300790 

QQ ( 5 * LT ) =QCALC ( AKK ( 5 ) ,XX,SX<4) ) 3 00800 

79 HL(5,LT)=CX(5)-SX(4)*QQ(5,LT) 300810 

85 CONTINUE 300820 

90 CONTINUE 300830 

95 RETURN 300840 

END 300850 


SUBROUTINE WHOUT(KDSTP) 

C BASIC OUTPUT AND INITIALIZATION 

COMMON /CHIN/ AK* AKK(ll), CNKH, LBEG, LLST • LL2 * INCtlO), 
1 KDSTP, KTQH * NSAV, NTOT » TIM(IO), TMPD 
COMMON /RUNWH/ AKH(20), AKAPT (20 ) * AREA ( 10 ) * DELT * 

1 HL (11*20)* HR (11*20)* KST(ll), KTRL(IO), NINT, NPTS, 

2 QQ (11*20) * SX<10), SY(10>, TMEI20) 

DIMENSION VEL ( 2 0 ) 

WR I TE ( 6 * 30 ) (TME(L) ,L = LBEG,LLST) 

WR I TE ( 6 ♦ 40 ) (AKAPT(L) , L=LBEG , LLST ) 

WR I TE ( 6 * 50 ) <AKH(L) ,L=LBEG,LLST) 

IL=NINT+1 
DO 28 1=1 tIL 

WR I TE ( 6 ,60 ) I» (QQ( I ,L ) »L=LBEG,LLST) 

DO 8 L = LBEG * LLST 
8 VEL (L) =QQ( I *L) /AREA! 1 ) 

WR I TE ( 6 * 70 ) (VEL(L) ,L=LBEG,LLST) 

I F ( 1-1)12.14,12 

12 WR I TE ( 6 » 8 0 ) ( HL ( I » L ) , L = LBEG , LLST ) 

14 IF< I-IL) 16,18.16 

16 WR I TE ( 6 ♦ 90 ) ( HR ( I , L ) , L = LBEG , LLST ) 

18 GO TO (22,28) .KDSTP 
22 LA=NSAV+1 

LB=LLST-NSAV 
DO 24 L=2»LA 
LB=LB+1 

QQ ( I . L ) =QQ ( I , LB ) 

HL ( I . L ) =HL ( I , LB ) 

HR ( I , L ) =HR ( I , LB ) 

24 CONTINUE 
28 CONTINUE 

30 FORMAT ( 10H1 TIM = 8G15.6) 

40 FORMAT ( 10HJ AKAPT = 8G15.6) 

50 FORMAT ( 10HJ AKH = 8G15.6) 

60 FORMAT ( 2HK ,I2,6H Q = 8G15.5) 

70 FORMAT ( 10HJ VEL = 8G15.6) 

80 FORMAT ( 10HJ HL = 8G15.6) 

90 FORMAT ( 10HJ HR = 8G15.6) 

RETURN 

END 
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FUNCTION QCALC (AKK,XX,YY) 500010 
WRK = 4.*ABS< XX ) / (AKK*AKK*YY*YY ) 5 00020 
WRK=SQRT( l.+WRK)-l. 500030 
WRK=AKK*AKK*YY*wRK/2. 500040 
QCALC = SIGN(WRK,XX ) 500050 
RETURN 500060 
END 500070 
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